Free Vibration Analysis of Curvilinearly Tapered Axially Functionally Graded Material Beams

1. Introduction

Incorporating variable cross-sections in structural elements facilitates achieving an optimal balance between weight and strength, especially in critical applications like helicopter blades, airplane propellers, turbine blades, and wind turbines. Given their significance, there has been extensive research on the vibrations of these structures, often modeled as beams with variable cross-sections [1,2,3].

The increasing use of functionally graded materials (FGMs) in beam-like structures underscores the need to account for material variations. FGMs enhance performance and provide a continuous stress distribution, unlike the discontinuities found in laminated and sandwich composites [4,5,6].

While initial research focused on material gradients in the thickness direction, there is significant interest in longitudinally graded materials, known as axially functionally graded material (AFGM) structures. Analyzing the free vibration of non-uniform AFGM beams involves recognizing that the coefficients in the equations of motion, derived from any beam theory, vary with changes in cross-sections and material properties along the beam's length. Consequently, solving these differential equations is not straightforward and often necessitates numerical procedures.

In the context of one-dimensional continuous models with inhomogeneous material parameters, Elishakoff and co-workers have presented exact solutions for the fundamental frequency of diverse AFG beams with different end supports [7]. Yuan et al. (2016) [8] applied the confluent hypergeometric function method to obtain exact solutions for Timoshenko beams with varying bending stiffness. Zhao et al. (2017) [9] employed the Chebyshev polynomials theory to analyze the free vibration of both AFGM Euler–Bernoulli and Timoshenko beams with non-uniform cross-sections. Xie et al. (2017) [10] exploited the spectral collocation method to examine the statics and free vibrations of Euler–Bernoulli beams with axially variable cross-sections and material parameters, while Chen (2021) [11] extended this technique to non-uniform AFG Timoshenko beams. Additionally, Zhang et al. (2019) [12] utilized a Jacobi polynomial-based approximation, and Cao et al. (2019) [13] used the asymptotic perturbation method for vibrational analysis of non-uniform AFGM beams.

Other researchers have creatively adapted methods known from studying homogeneous beams with variable cross-sections to address solutions for non-uniform AFGM beams. For instance, Ghazaryan et al. (2018) [14] employed the differential transform method to study the free vibration characteristics of AFGM Euler-Bernoulli beams, while Rajasekaran and Tochaei (2014) [15] used both the differential transform element and differential quadrature element methods to study AFGM Timoshenko beams, both with varying cross-sectional profiles and material properties along the beam axis. Keshmiri et al. (2018) [16] determined natural frequencies of cantilever beams with different variable cross-sections and material gradients using the Adomian decomposition method (ADM), and Lin et al. (2022) [17] investigated the vibration of rotating non-uniform AFG Euler–Bernoulli beams using the Laplace ADM. In turn, Wang et al. (2022) [18] combined ADM with an iterative process for reliability analysis of composite beams with varying cross-section rigidities and mass distributions along their length. Mahmoud (2019) [19] presented a general solution for the free transverse vibration of cantilever non-uniform AFGM Euler-Bernoulli beams, loaded at the tips with point masses, using the Myklestad method. Chen et al. (2021) [20] employed the variational iteration method to determine the modal characteristics of tapered AFGM Euler–Bernoulli and Timoshenko beams. Liu et al. (2022) [21] proposed a closed-form dynamic stiffness formulation for free vibration analysis of tapered and/or FGM Euler–Bernoulli beams. Adelkhani and Ghanbari (2022) [22] analyzed the vibrations of tapered AFGM beams with nonlinear profiles using the point collocation method.

To investigate the dynamics of AFGM tapered beams under large deflection scenarios, Kumar et al. (2015) [23] employed the Rayleigh-Ritz method with start functions derived from nonlinear static analysis within the Euler-Bernoulli beam theory with von Kármán geometric nonlinearity, while Ghayesh (2018) [24] utilized the Galerkin method and the third-order shear deformation beam theory. Soltani and Asgarian (2019) [25] developed a hybrid approach combining the power series expansions and the Rayleigh-Ritz method for stability and free vibration analyses of non-uniform AFGM beams resting on Winkler-Pasternak elastic foundation. Singh and Sharma (2022) [26] carried out a vibration analysis of an AFGM non-prismatic Timoshenko beam under axial thermal variation in a humid environment using the harmonic differential quadrature method.

Finite element method-based approaches, tailored for precise and efficient study of this task, have also been proposed. Özdemir (2022) [27] studied free vibration and buckling characteristics of AFGM tapered rotating Euler–Bernoulli and Timoshenko beams with different end conditions using a two-noded beam element. Bazoune (2024) [28] proposed a Fourier-p element model for accurate predictions of Timoshenko-Ehrenfest beam frequencies with high precision compared to other existing models. Chen et al. (2019) [29] presented isogeometric analysis (IGA) in conjunction with elasticity theory for the three-dimensional vibration problem of AFGM beams with variable thickness. Murillo et al. [30] used the seven-parameter spectral finite element formulation to perform the analysis of FG shells of either uniform or non-uniform thickness. Burlayenko et al. (2024) [31] conducted simulations on the free vibrations of AFGM beams with non-uniform cross-sections, utilizing both one-dimensional and three-dimensional AFG finite element models developed in ABAQUS via user-defined subroutines [32,33,34].

Most of the examined non-uniform beams had rectangular cross-sections with linear taper. There were also studies on beams with parabolic and exponential thickness and constant or linear varying width, but less common. In some cases, vibrations of beams were investigated in scenarios such as truncated beams, beams with one sharp end, and beams with both ends sharp, as discussed in [1,3]. However, research on AFGM beams having arbitrary cross-section variations remains limited. Recently, Lee and Lee (2022) [35] studied the coupled flexural-torsional free vibration of circular horizontally curved beams with rectangular and elliptical cross-sections made of AFGMs with quadratic functions of Young’s modulus and the mass density, using the trial eigenvalue method and numerical integration within the Timoshenko and Saint-Venant beam theories. Rezaiee-Pajand et al. (2022) [36] proposed a closed-form solution for the lateral-torsional buckling moment of a bidirectional exponentially functionally graded monosymmetric C-shaped Euler–Bernoulli beam. Liu et al. (2024) [37] presented a nonlinear model to obtain the mode shape and frequency of FGM Euler-Bernoulli thin-walled beams with varying cross-sections, employing the perturbation approach and the Galerkin method. Other researchers have explored the concept of transforming the equations of motion of non-uniform FGM beams into those of equivalent uniform beams. For instance, Chen et al. (2017) [38] applied this idea to a specific class of non-uniform FGM beams. Additionally, Martin and Salehian (2020) [39] developed a method based on metric minimization, modal participation factor and the proper orthogonal decomposition to accurately approximate a Euler-Bernoulli motion equation with spatially varying coefficients with an equivalent constant coefficient model.

The literature search indicates that studying the dynamics of tapered AFGM beams involves complexities. While various techniques have been reported, closed-form solutions are limited to specific cross-sections, material variations and defined boundary conditions. Numerical methods like FEM are proficient but make it challenging to derive general conclusions about system behavior [40]. In this respect, accurate semi-analytical approaches can be a good alternative. Despite extensive research, the vibration characteristics of AFGM beams with arbitrary continuously variable cross-sections are not fully addressed. This study aims to use the differential transform method (DTM) to obtain accurate semi-analytical solutions for the free vibration of curvilinearly tapered AFGM beams under various boundary conditions, determining the effects of variable cross-sections on natural frequencies and mode shapes. This is an innovation of the present contribution.

2. AFGM Beams with Variable Cross-Section

2.1. Geometry of the Beams

Focusing on the closed-form solution of the free vibration problem, we aim to describe various non-uniform beam geometries in a unified analytical manner. One effective way to achieve this is by using rational Bézier curves [41].

We consider planar beams of length of $L$ that are symmetrical about their midlines, with the left ends positioned at the origin of the Cartesian coordinate system as shown in Figure 1. The top curve of these beams is a conic arc, represented by an algebraic quadratic equation in ${\mathbb{R}}^2$ as follows:

$$A{x}^2+Bxy+C{x}^2+Dx+Ey+F=0.$$

(1)

The coefficients of (1) can be derived from the quadratic rational Bézier curve, which involves the control points $P_0(x_0;y_0)$, $P_1(x_1;y_1)$, $P_2(x_2;y_2)$ and the weight $w_1$, by utilizing the Gröbner basis as presented in [42]:

$$\begin{array}{l}A=y_0^2-4w_1^2{y}_0{y}_1+4w_1^2{y}_1^2-2y_0{y}_2+4w_1^2{y}_0{y}_2-4w_1^2{y}_1{y}_2+y_2^2,\\ B=-2x_0{y}_0+4w_1^2{x}_1{y}_0+2x_2{y}_0-4w_1^2{x}_2{y}_0+4w_1^2{x}_0{y}_1-8w_1^2{x}_1{y}_1+4w_1^2{x}_2{y}_1+2x_0{y}_2-4w_1^2{x}_0{y}_2+4w_1^2{x}_1{y}_2-2x_2{y}_2,\\ C=x_0^2-4w_1^2{x}_0{x}_1+4w_1^2{x}_1^2-2x_0{x}_2+4w_1^2{x}_0{x}_2-4w_1^2{x}_1{x}_2+x_2^2,\\ D=-2x_2{y}_0^2+4w_1^2{x}_1{y}_0{y}_1+4w_1^2{x}_2{y}_0{y}_1-4w_1^2{x}_0{y}_1^2-4w_1^2{x}_2{y}_1^2+2x_0{y}_0{y}_2-8w_1^2{x}_1{y}_0{y}_2+2x_2{y}_0{y}_2+4w_1^2{x}_0{y}_1{y}_2+4w_1^2{x}_1{y}_1{y}_2-2x_0{y}_2^2,\\ E=-4w_1^2{x}_1^2{y}_0+2x_0{x}_2{y}_0+4w_1^2{x}_1{x}_2{y}_0-2x_2^2{y}_0+4w_1^2{x}_0{x}_1{y}_1-8w_1^2{x}_0{x}_2{y}_1+4w_1^2{x}_1{x}_2{y}_1-2x_0^2{y}_2+4w_1^2{x}_0{x}_1{y}_2-4w_1^2{x}_1^2{y}_2+2x_0{x}_2{y}_2,\\ F=x_2^2{y}_0^2-4w_1^2{x}_1{x}_2{y}_0{y}_1+4w_1^2{x}_0{x}_2{y}_1^2+4w_1^2{x}_1^2{y}_0{y}_2-2x_0{x}_2{y}_0{y}_2-4w_1^2{x}_0{x}_1{y}_1{y}_2+x_0^2{y}_2^2.\end{array}$$

(2)

For convenience, let $l_f$ and $l_e$ denote the length of lines ${\mathit{P}}_1{\mathit{P}}_m$ and ${\mathit{P}}_s{\mathit{P}}_m$, respectively, where ${\mathit{P}}_S=\left(1-s\right){\mathit{P}}_m+s{\mathit{P}}_1$ is the shoulder point, ${\mathit{P}}_m$ is the midpoint of line ${\mathit{P}}_0{\mathit{P}}_2$ and $s={\displaystyle \frac{l_e}{l_f}}$ is a shape parameter related to $w_1={\displaystyle \frac{s}{1-s}}$. By specifying the coordinates of the points ${\mathit{P}}_0$, ${\mathit{P}}_1$, ${\mathit{P}}_2$ and the parameter $s$, the implicit equation of a conic curve (1) can be constructed. The curve will have a convex upward shape if the vertex ${\mathit{P}}_1$ is above the line ${\mathit{P}}_0{\mathit{P}}_2$, and a concave downward shape if ${\mathit{P}}_1$ below the line ${\mathit{P}}_0{\mathit{P}}_2$ as illustrated in Figure 1a and b, respectively.

2.2. Graded Materials of the Beams

We suppose that the materials of the beams with variable thickness are homogeneous in the thickness direction but graded along the axial direction. Consequently, the material properties, such as mass density and Young’s modulus, vary with the position along the $x$-axis. Specifically, we analyze the following types of material distributions:

Type 1. The effective mass density and Young’s modulus are power-law functions of $x$ written as

$$\rho \left(x\right)={\rho }_0{\sum _{i=0}^m{a}_i\left({\displaystyle \frac{x}{L}}\right)}^i,\hspace{1em} E\left(x\right)=E_0{\sum _{i=0}^n{b}_i\left({\displaystyle \frac{x}{L}}\right)}^i,$$

(2)

where ${\rho }_0$ and $E_0$ are the mass density and Young’s modulus at $x=0$, respectively; $m$ and $n$ are power coefficients that account for arbitrary variations of the material parameters along the beam’s length; $a_i$ and $b_i$ constant coefficients.

Type 2. The effective density and Young’s modulus follow exponential law functions of $x$, given by

$$\rho \left(x\right)={\rho }_0{V}_e, E\left(x\right)=E_0{V}_e,$$

(3)

where the exponential function $V_e=e^{{\gamma }_e\left(1-a\left(1-{\displaystyle \frac{x}{L}}\right)\right)}$ is determined by the parameter $a$ and gradation index ${\gamma }_e$, [29].

Type 3. The beams with variable thickness are composed of a metal-ceramic mixture. The effective material properties are calculated based on the rule of mixtures in the form:

$$\rho \left(x\right)={\rho }_c+\left({\rho }_m-{\rho }_c\right)V_c, E\left(x\right)=E_c+\left(E_m-E_c\right)V_c,$$

(4)

where the subscripts ‘$m$’ and ‘$c$’ refer to the properties of metal and ceramic constituents, respectively; $V_c$ denotes the ceramic volume fraction, expressed as $V_c={\left(1-b\left(1-{\displaystyle \frac{x}{L}}\right)\right)}^{{\gamma }_c}$ with the parameter $b$ and gradation index ${\gamma }_c$, [29].

Poisson’s ratio, $\nu$ is assumed to be constant throughout this study.

2.3. Governing Equations of the Beams

We consider the transverse vibrations of an AFGM beam in the $xz$-plane using the Timoshenko beam theory. The equations of motion governing such beams with generally variable cross-sections, in terms of the transverse displacement $w(x,t)$ of the midline and the cross-sectional rotation angle $\phi (x,t)$ at time $t$ are written as follows:

$$\begin{gathered} -{\displaystyle \frac{\partial }{\partial x}}\left(D_{11}\left(x\right){\displaystyle \frac{\partial \phi }{\partial x}}\right)-A_{55}\left(x\right)\left({\displaystyle \frac{\partial w}{\partial x}}-\phi \right)+I_2\left(x\right){\displaystyle \frac{{\partial }^2\phi }{\partial t^2}}=0, \\ -{\displaystyle \frac{\partial }{\partial x}}\left(A_{55}\left(x\right)\left({\displaystyle \frac{\partial w}{\partial x}}-\phi \right)\right)+I_0\left(x\right){\displaystyle \frac{{\partial }^{2}w}{\partial t^2}}=0. \end{gathered}$$

(5)

The rigidities and the inertia coefficients are expressed as integrals over the cross-sectional domain in the form:

$$\begin{gathered} D_{11}\left(x\right)=\underset{A\left(x\right)}{\overset{ }{\int }}{z}^{2}E\left(x\right)dA, A_{55}\left(x\right)=\kappa \underset{A\left(x\right)}{\overset{ }{\int }}G\left(x\right)dA, \\ \left(I_0\left(x\right),I_2\left(x\right)\right)=\underset{A\left(x\right)}{\overset{ }{\int }}\left(1,z^2\right)\rho \left(x\right)dA, \end{gathered}$$

(6)

where $\rho \left(x\right)$ is the mass density, $G\left(x\right)={\displaystyle \frac{E\left(x\right)}{2(1+\nu )}}$ is the shear modulus and $\kappa$ is the shear correction factor. Throughout this study, the shear factor is accepted as 5/6 unless specified otherwise.

Assuming simple harmonic oscillations, the deflection and the rotation angle can be expressed as

$$w\left(x,t\right)=\overline{w}\left(x\right)e^{i\omega t}, \phi \left(x,t\right)=\overline{\phi }\left(x\right)e^{i\omega t}$$

(7)

where $\overline{w}\left(x\right)$ and $\overline{\phi }\left(x\right)$ are the amplitudes of the displacement and rotation angle, respectively, $\omega$ is the angular frequency; and $i^2=-1$ represents the imaginary unit.

Substituting (7) into (5) yields a system of coupled ordinary differential equations that describe the free transverse vibrations of a Timoshenko beam with a non-uniform cross-section as

$$\begin{gathered} {\displaystyle \frac{d}{dx}}\left(D_{11}\left(x\right){\displaystyle \frac{d\overline{\phi }}{dx}}\right)+A_{55}\left(x\right)\left({\displaystyle \frac{d\overline{w}}{dx}}-\overline{\phi }\right)+I_2\left(x\right){\omega }^2\overline{\phi }=0, \\ {\displaystyle \frac{d}{dx}}\left(A_{55}\left(x\right)\left({\displaystyle \frac{d\overline{w}}{dx}}-\overline{\phi }\right)\right)+I_0\left(x\right){\omega }^2\overline{w}=0. \end{gathered}$$

(8)

Herewith, the rigidities and inertia coefficients (6) can be conveniently presented by the expressions:

$$\begin{gathered} D_{11}\left(x\right)=E\left(x\right)I\left(x\right), A_{55}\left(x\right)=\kappa G\left(x\right)A\left(x\right) \\ {I}_0\left(x\right)=\rho \left(x\right)A\left(x\right), I_2\left(x\right)=\rho \left(x\right)J\left(x\right), \end{gathered}$$

(9)

where$A\left(x\right)$ and $J\left(x\right)$ are the area and second inertia moment of the cross-section, respectively.

3. Solution Methodology

3.1. Differential Transform Method

Due to the variable coefficients, reflecting the dependence of material properties and cross-sectional geometry on the $x$-coordinate, the coupled nature of the differential equations in (8) presents mathematical challenges for finding a closed-form solution. To address this, we employ the differential transform method (DTM) for its versatility and simplicity. This method converts the differential equations into a set of recurrent algebraic equations, enabling the solution to be represented as an infinite power series. The fundamental aspects of the DTM can be found in G.E. Pukhov’s original publications [43].

We reformulate system (8) a more suitable form for applying the DTM as follows:

$$\begin{gathered} {\displaystyle \frac{d^2\overline{\phi }}{d{x}^2}}+{\overline{D}}_1\left(x\right){\displaystyle \frac{d\overline{\phi }}{dx}}+{\overline{D}}_2\left(x\right){\displaystyle \frac{d\overline{w}}{dx}}-{\overline{D}}_2\left(x\right)\overline{\phi }+{\omega }^2{\overline{D}}_3\left(x\right)\overline{\phi }=0, \\ {\displaystyle \frac{d^2\overline{w}}{d{x}^2}}+{\overline{G}}_1\left(x\right){\displaystyle \frac{d\overline{w}}{dx}}-{\displaystyle \frac{d\overline{\phi }}{dx}}-{\overline{G}}_1\left(x\right)\overline{\phi }+{\omega }^2{\overline{G}}_2\left(x\right)\overline{w}=0, \end{gathered}$$

(10)

where the variable coefficients are defined as ${\overline{D}}_1\left(x\right)={\displaystyle \frac{D_{11}^{\prime }\left(x\right)}{D_{11}\left(x\right)}}$, ${\overline{D}}_2\left(x\right)={\displaystyle \frac{A_{55}\left(x\right)}{D_{11}\left(x\right)}}$, ${\overline{D}}_3\left(x\right)={\displaystyle \frac{I_2\left(x\right)}{D_{11}\left(x\right)}}$, ${\overline{G}}_1\left(x\right)={\displaystyle \frac{A_{55}^{\prime }\left(x\right)}{A_{55}\left(x\right)}}$, ${\overline{G}}_2\left(x\right)={\displaystyle \frac{I_0\left(x\right)}{A_{55}\left(x\right)}}$. Here, and in what follows, the prime denotes differentiation with respect to the $x$-coordinate.

Hence, the solution of (10) at a specific point $x_0$ within the interval $0\le x\le L$ in the domain of differential transformation is expressed by the recurrent algebraic equations:

$$\begin{array}{l}\mathsf{\Phi }\left(k+2\right)=k_{*}\{-\sum _{p=0}^k\left(p+1\right)\mathsf{\Phi }\left(p+1\right)D_1\left(k-p\right)+\sum _{p=0}^k\mathsf{\Phi }\left(p\right)D_2\left(k-p\right)\\ -\sum _{p=0}^k\left(p+1\right)W\left(p+1\right)D_2\left(k-p\right)-{\omega }^2\sum _{p=0}^k\mathsf{\Phi }\left(p\right)D_3\left(k-p\right)\}\\ W\left(k+2\right)=k_{*}\{-\sum _{p=0}^k\left(p+1\right)W\left(p+1\right)G_1\left(k-p\right)+\sum _{p=0}^k\mathsf{\Phi }\left(p\right)G_1\left(k-p\right)\\ +\sum _{p=0}^k\left(p+1\right)\mathsf{\Phi }\left(p+1\right)L\left(k-p\right)-{\omega }^2\sum _{p=0}^{k}W\left(p\right)G_2\left(k-p\right)\}\end{array}$$

(11)

where $k_{*}={\displaystyle \frac{k!}{\left(k+2\right)!}}$, $L\left(k\right)$ is a unity image, and $W\left(k\right)$, $\mathsf{\Phi }\left(k\right)$, $D_1\left(k\right)$, $D_2\left(k\right)$, $D_3\left(k\right)$, $G_1\left(k\right)$, and $G_2\left(k\right)$ are the images of the unknown amplitudes $\overline{w}\left(x\right)$ and $\overline{\phi }\left(x\right)$ and the corresponding functional coefficients in (10), respectively.

By incrementing the index $k$ sequentially, we can generate all images from (11) except for $W\left(0\right)$, $W\left(1\right)$, $\mathsf{\Phi }\left(0\right)$, and $\mathsf{\Phi }\left(1\right)$. Consequently, the system (11) can be simplified to the form:

$$\begin{gathered} \mathsf{\Phi }(k+2)=B_{\phi }\left(k\right)W\left(0\right)+C_{\phi }\left(k\right)\mathsf{\Phi }\left(0\right)+G_{\phi }\left(k\right)W\left(1\right)+H_{\phi }\left(k\right)\mathsf{\Phi }\left(1\right), \\ W(k+2)=B_w\left(k\right)W\left(0\right)+C_w\left(k\right)\mathsf{\Phi }\left(0\right)+G_w\left(k\right)W\left(1\right)+H_w\left(k\right)\mathsf{\Phi }\left(1\right) \end{gathered}$$

(12)

where the explicit forms of the recurrent coefficients in (12) are presented in [44] (see Appendix B).

Given that the images $W\left(k\right)$ and $\mathsf{\Phi }\left(k\right)$ at $k\ge 2$ depend on $W\left(0\right)$, $W\left(1\right)$, $\mathsf{\Phi }\left(0\right)$ and $\mathsf{\Phi }\left(1\right)$, the original functions $\overline{w}\left(x\right)$ and $\overline{\phi }\left(x\right)$ can be reconstructed for a selected number of images $N$ as follows:

$$\begin{gathered} \overline{w}\left(x,\omega \right)=\left\{1+\sum _{p=0}^{N-2}{B}_w\left(p\right){\left(x-x_0\right)}^{p+2}\right\}W\left(0\right)+\left\{\sum _{p=0}^{N-2}{C}_w\left(p\right){\left(x-x_0\right)}^{p+2}\right\}\mathsf{\Phi }\left(0\right) \\ +\left\{(x-x_0)+\sum _{p=0}^{N-2}{G}_w\left(p\right){\left(x-x_0\right)}^{p+2}\right\}W\left(1\right)+\left\{\sum _{p=0}^{N-2}{H}_w\left(p\right){\left(x-x_0\right)}^{p+2}\right\}\mathsf{\Phi }\left(1\right), \\ \overline{\phi }\left(x,\omega \right)=\left\{\sum _{p=0}^{N-2}{B}_{\phi }\left(p\right){\left(x-x_0\right)}^{p+2}\right\}W\left(0\right)+\left\{1+\sum _{p=0}^{N-2}{C}_{\phi }\left(p\right){\left(x-x_0\right)}^{p+2}\right\}\mathsf{\Phi }\left(0\right) \\ +\left\{\sum _{p=0}^{N-2}{G}_{\phi }\left(p\right){\left(x-x_0\right)}^{p+2}\right\}W\left(1\right)+\left\{(x-x_0)+\sum _{p=0}^{N-2}{H}_{\phi }\left(p\right){\left(x-x_0\right)}^{p+2}\right\}\mathsf{\Phi }\left(1\right). \end{gathered}$$

(13)

Taking into account the actual boundary conditions applied to the beam ends, which involve the transverse displacement $\overline{w}$, rotation angle $\overline{\phi }$, bending moment $\overline{M}$, and shear force $\overline{Q}$, such that

$$\overline{M}=D_{11}\left(x\right){\displaystyle \frac{d\overline{\phi }}{dx}}, \overline{Q}=A_{55}\left(x\right)\left({\displaystyle \frac{d\overline{\overline{w}}}{dx}}-\overline{\phi }\right)$$

we can formulate the eigenvalue problem for each case of constraints in the form:

$$\left[\begin{array}{cccc}{A}_{11}\left(\omega \right)& A_{12}\left(\omega \right)& A_{13}\left(\omega \right)& A_{14}\left(\omega \right)\\ A_{21}\left(\omega \right)& A_{22}\left(\omega \right)& A_{23}\left(\omega \right)& A_{24}(\omega \\ A_{31}\left(\omega \right)& A_{32}\left(\omega \right)& A_{33}\left(\omega \right)& A_{34}(\omega \\ A_{41}\left(\omega \right)& A_{42}\left(\omega \right)& A_{43}\left(\omega \right)& A_{44}(\omega \end{array}\right]\left(\begin{array}{c}W\left(0\right)\\ \mathsf{\Phi }\left(0\right)\\ W\left(1\right)\\ \mathsf{\Phi }\left(1\right)\end{array}\right)=0.$$

(14)

Here, the functions $A_{ij}\left(\omega \right)$ are appropriate polynomials of the unknown angular frequency $\omega$, which involve the recurrent coefficients and the power terms of $(x – x_0)$.

Finally, the eigenvalues are computed by solving the root-finding problem resulting from the equation:

$$\det \left|A_{ij}\left(\omega \right)\right|=0.$$

(15)

A computational program has been developed in the Matlab environment to apply the DTM for analyzing the free vibration of curvilinear tapered inhomogeneous Timoshenko beams. The eigenvalue problem is addressed using standard algorithms provided in the Matlab package [45].

3.2. Finite Element Modeling

To ensure the accuracy and validation of the proposed semi-analytical methodology for analyzing the free vibration parameters of curvilinearly tapered AFGM beams, a finite element model has been developed and utilized for comparison purposes. It is well-known that the optimal modeling strategy involves using graded finite elements (GFEs), which allow for the variation of material properties at the element level. Additionally, advanced three-dimensional models are essential to precisely capture the non-uniform geometry of axially inhomogeneous beams [31].

The free vibration analysis of AFGM beams is conducted using ABAQUS [46], incorporating a user-defined material subroutine (UMAT) to implement material gradients into a convenient three-dimensional element, C3D20, as detailed in [31,33] and illustrated in Figure 2b. The material gradient is assigned within the element volume during the computation of the element stiffness matrix, ${\mathit{K}}^e$ using Gauss quadrature, as follows:

$${\mathit{K}}^e={\sum }_{p=1}^{n_p}{\sum }_{q=1}^{n_q}{\sum }_{r=1}^{n_r}{\left({\mathit{B}}^T\mathit{D}\left(x\right)\mathit{B}\left|\mathit{J}\right|\right)}_{p,q,r}{w}_p{w}_q{w}_r,$$

(16)

where the indices $p$, $q$, and $r$ refer to the integration points over the element volume, $\left|\mathit{J}\right|$ is the determinant of the geometric Jacobian matrix, $w_p$, $w_q$, and $w_r$ are the Gaussian weights; $\mathit{B}$ stands for a matrix of gradients of the shape functions for the 3-D element and the material matrix $\mathit{D}\left(x\right)$ of the 3-D GFE relevant to this study can be found in [33].

ABAQUS currently does not support calculating the mass matrix based on a user-defined mass density in frequency analysis [46]. The elements available in ABAQUS use a constant mass density, i.e.

$${\mathit{M}}^e={\int }_{V^e}{\rho }_0{\mathit{N}}^T·\mathit{N}d{V}^e,$$

(17)

where $\mathit{N}$ is an element interpolation matrix.

To address this limitation, we assigned a constant effective mass density within the GFE by using an average value of the actual mass density distribution over the beam's volume in the form:

$${\rho }_{avr}={\displaystyle \frac{1}{V}}{\int }_V\rho \left(x\right)dx$$

(18)

In the case of the free vibration analysis without damping, the global discrete system of governing equations of motion is reduced to the problem of eigenvalues and associated eigenvectors, which is written as follows

$$\left(\mathit{K}-{\omega }^2\mathit{M}\right)\mathsf{{\rm X}}=0$$

(19)

where $\mathit{M}$, $\mathit{K}$ are global mass and stiffness matrices, obtained by assembling the formulated element mass and stiffness matrices; $\omega$ is the undamped circular frequency, $\mathsf{{\rm X}}$ is the vector of nodal displacements associated with oscillations at each $\omega$. The free vibration analysis is carried out using linear perturbation load step, obtaining the eigenvalues by means of Lanczos or subspace iteration methods [46]. Figure 2a illustrates a 3-D finite element beam model created using the ABAQUS pre-processor, utilizing the 'loft' feature. To provide accurate results, the element aspect ratio in the mesh was consistently maintained close to one throughout all subsequent calculations.

4. Results

We will now examine beams with various curvilinear tapers in rectangular cross-sections, incorporating three distinct material variations along the beam axis. This analysis aims to demonstrate the effectiveness of the proposed computational technique by comparing the results with those obtained from other methods. First, we validate the accuracy of the DTM approach by benchmarking our findings against known solutions for uniform homogeneous and inhomogeneous Timoshenko beams under various boundary conditions. Then, we present a selection of examples for curvilinear taper beams, including both known and previously unexplored cases, to further illustrate the method's applicability and provide new insights that may be valuable to other researchers in this field.

4.1. Homogeneous Timoshenko Beams with Uniform Cross-Sections

The accuracy and effectiveness of the proposed approach are first benchmarked by computing the frequencies of homogeneous beams with uniform cross-sections.

We consider geometrically uniform ($w_1=0$) homogeneous ($m=n=0$ with unit coefficients in (2) ) beams with clamped-clamped (C-C), clamped-free (C-F), and simply supported (S-S) boundary conditions. The material and geometric parameters of the beams are listed in Table 1.

The first seven natural frequencies of the S-S, C-C, and C-F beams compared with those computed using the dynamic stiffness matrix method (DSM) in conjunction with the Timoshenko beam theory in [47] are presented in Table 2.

As shown in Table 2, there is excellent agreement between the solutions obtained using the DTM and DSM methods for all vibrational modes.

4.2. Inhomogeneous Timoshenko beams with uniform cross-sections

Next, the ability of the proposed approach to accurately predict the frequencies of beams with uniform cross-sections made of inhomogeneous materials is evaluated.

We assume that the beam’s mass density and modulus of elasticity follow power-law distributions of Type 1 as described in (2). Since natural frequencies for Timoshenko beams under these conditions are not documented in the literature, we will compare our results with those obtained using the Euler-Bernoulli beam theory in [14] for AFGM beams with a slenderness ratio of $L/h=50$. The geometrical and material parameters are consistent with those detailed in [14] (see Tables 4–6 and Appendix B).

The comparison of the first five dimensionless frequencies $\overline{\omega }=\omega L^2\sqrt{{\displaystyle \frac{{\rho }_{0}A}{E_{0}J}}}$ for the cantilever AFGM beam is presented in Table 3. Furthermore, to extend this study to beams, where shear deformation plays a significant role, we also provide in Table 3 the natural frequencies computed for beams with a slenderness ratio of 10.

Upon reviewing Table 3, it becomes evident that both Euler-Bernoulli and Timoshenko beams with slenderness ratio of 50 demonstrate excellent agreement for the fundamental frequency and show slight deviations for higher frequencies. This is because the Euler-Bernoulli beam theory overestimates especially higher frequencies due to its neglect of shear effect. Specifically, when the slenderness ratio is 10, the fundamental frequencies of the beams show only satisfactory matching, while significant discrepancies are observed in higher-order frequencies, consistent with theoretical expectations.

In the second example, we evaluate the accuracy of predicted natural frequencies for exponential AFGM beams. The current results are compared with those obtained in [11] using the Chebyshev collocation method. The beams are composed of zirconia and aluminum with axial gradation according to law (4), where the ceramic volume fraction is an exponential function: $V_c={\displaystyle \frac{e^{\alpha x/L}-1}{e^{\alpha }-1}}$. The material and geometrical properties of the beams are described in [11] and detailed in Table 4 as follows:

The computed first four dimensionless natural frequencies $\overline{\omega }=\omega L^2\sqrt{{\displaystyle \frac{{\rho }_{Al}A}{E_{Al}J}}}$ of the exponential-law AFGM beams under S-S and C-C boundary conditions for various exponents $\alpha$ compared with the results in [11] are given in Table 5.

From Table 5, it is obvious that the present results for all first four frequencies closely align with those obtained in the referenced study.

Finally, uniform aluminum-zirconia Timoshenko beams with axial variation of the material characteristics according to the power-law function $V_c={\left({\displaystyle \frac{x}{L}}\right)}^{\beta }$ of the ceramic volume fraction in the rule of mixtures (4) are tested. The natural frequencies of these beams are benchmarked against the findings in [15]. The parameters of the AFGM beams are the same as those in [15] and are listed in Table 4.

The first four dimensionless natural frequencies $\overline{\omega }=\omega L^2\sqrt{{\displaystyle \frac{{\rho }_{Zr}A}{E_{Zr}J}}}$ of the AFGM beams under C-F, S-S and C-C boundary conditions for various indices $\beta$, compared with the solutions obtained using the differential quadrature element method in [15], are shown in Table 6.

Inspecting Table 6 reveals excellent agreement between the first four frequencies obtained in this study and those reported in the referenced research across all boundary conditions and power indices.

It should be noted that throughout this series of examples, the proposed approach demonstrates high accuracy, efficiency, and rapid convergence. Only thirty terms in the series expansions (13) are sufficient to yield accurate results for the first several natural frequencies.

4.3. Homogeneous Timoshenko Beams with Curvilinearly Tapered Cross-Sections

In order to demonstrate the applicability and effectiveness of the DTM in solving the free vibration problem for axially curvilinear cross-sectional beams, we now consider various specific cases. Whenever possible, we compare the results obtained using DTM with those already documented in the literature. In cases, where comparative data is not available, we perform the frequency analysis using ABAQUS software.

The first problem addressed is a homogeneous beam with a double-parabolic taper, truncated at the left end. The shape of the parabolic taper beam is analytically defined by the weight $w_1=1$, the parameter $l_e$ equals to the truncation factor ${\delta }_0$ and the tringle ${\mathit{P}}_0{\mathit{P}}_1{\mathit{P}}_2$ specified as shown in Figure 3a.

Due to the lack of data for such Timoshenko beams in the published literature, the results are compared with those obtained using the Rayleigh-Ritz method within the Euler-Bernoulli beam theory in [2]. Table 7 presents the first three dimensionless frequencies $\overline{\omega }=\omega {\left(L-{\delta }_0\right)}^2\sqrt{{\displaystyle \frac{\rho A_0}{E{J}_0}}}$ for both the reference Euler-Bernoulli and present Timoshenko beams with a slenderness ratio $L/h_0=50$ under three different boundary conditions and different truncation factors.

One can observe from Table 7 that the eigenfrequencies closely match for both types of beams when the truncation factor is small. However, as the truncation factor increases, discrepancies between them become more pronounced, especially for the third mode. This difference is anticipated because shorter beams necessitate accounting for shear deformation, which the simpler Euler-Bernoulli beam model does not address.

The second benchmark example involves a beam with a parabolic-tapered height and constant width, truncated at the left end and tapered to a point at the right end (${\delta }_1=0$), as illustrated in Figure 3b. A detailed specification of this beam configuration is provided in [3]. The analytical presentation of this shape of the parabolic beam with Bézier curves is discussed in Section 2.1.

Once again, due to limited availability of data for Timoshenko beams in the existing literature, comparisons are made with results obtained using hypergeometric functions within the framework of the Euler-Bernoulli beam theory, as described in [3]. Table 8 displays the first four dimensionless frequencies $\overline{\omega }=\omega {\left(L/2\right)}^2\sqrt{{\displaystyle \frac{\rho A_0}{E{J}_0}}}$ of cantilever Euler-Bernoulli and Timoshenko beams with a slenderness ratio $L/h_0=50$ for various truncation factors.

From Table 8, one can see that the compared natural frequencies exhibit close agreement. Similar to the previous example, as the beam length decreases, discrepancies between the present and benchmarked results increase, especially for higher modes, due to the greater influence of shear deformation, as discussed earlier.

Therefore, these examples thoroughly confirm the reliability of the proposed computational method for analyzing curvilinear tapered beams.

4.4. AFGM Timoshenko Beams with Curvilinearly Tapered Cross-Sections

This subsection presents the performance of the DTM for frequency analysis of inhomogeneous curvilinearly tapered beams of convex-up and concave-down shapes, referred to as plump and slender models in [29], illustrated in Figure 1a,b, respectively. The material properties of the beams vary axially according to the laws specified in (3) for Type 2 and (4) for Type 3 distributions. Detailed geometric and material properties of the beams can be found in [29] and are summarized in Table 9.

First, we analyze beams with varying thickness and constant width for both Type 2 and Type 3 material laws. The exponential law (3) with parameter $a=1$ and gradient indices ${\gamma }_e=0, 1, 5$ is applied to beams made of aluminum, whereas the power-law (4) with parameter $b=1$ and gradient indices ${\gamma }_c=0, 1, 5$ is utilized for beams composed of a combination of steel (SUS304) and ceramic (Si₃N₄) phases, as done in [29].

Table 10 presents the first two frequencies associated with the bending modes of the curvilinear taper beams with C-C boundary conditions. These frequencies are compared with those computed using NURBS in the 3-D finite element analysis conducted in [29].

The current results align satisfactory with the reference data, demonstrating differences of less than 0.6% for the fundamental frequency and less than 9% for the second frequency in the plump beams. For slender beams, the discrepancies are approximately 3% for both frequencies. These deviations arise from differences in modeling approaches. The current model uses a line beam structure with a first-order shear theory, whereas the IGA analysis in [29] employs a more detailed 3-D beam modeling approach.

Next, we analyze the parabolic convex beams depicted in Figure 3b. These beams feature a parabolic taper in height and maintain a constant width. They are truncated at the left end with a fixed factor ${\delta }_0=0.5$, while the right end is characterized by a variable truncation factor ${\delta }_1$. The beam ends are subjected to C-C and C-F constraints. The material constants, given in Table 1, vary along the beam length according to the Type 1 material law with a linear distribution for Young's modulus ($n=1$) and a parabolic distribution for mass density ($m=2$) as described in (2). Detailed specifications of these beam configurations are provided in [23].

While benchmark results for the fundamental frequencies of these beams are available in [23], higher-order frequencies are absent in the existing literature. Therefore, the results obtained using the DTM are validated against finite element analysis for these beams, as discussed in Section 3.2.

The comparison of first five dimensionless bending frequencies $\overline{\omega }=\omega L^2\sqrt{{\displaystyle \frac{{\rho }_0{A}_0}{E_0{J}_0}}}$, depending on the truncation factor ${\delta }_1$ for the two types of boundary conditions is shown in Table 11. It is observed that the differences between dimensionless fundamental frequencies obtained from the present analysis, published results and FEM predictions are quite small. However, the higher frequencies show more significant deviations between the DTM and FEM solutions. Likely, these discrepancies between the two sets of data are attributed to methodological differences in formulation and solution approach. Overall, these findings confirm the reliability of the DTM approach within acceptable limits.

As mode shapes play a crucial role in evaluating the accuracy of structural dynamics modeling, the first five bending mode shapes associated with the appropriate frequencies of AFGM beams with ${\delta }_1=0.2$ for C-C and C-F boundary conditions, as listed in Table 11 and obtained using the DTM, are plotted in Figure 4a,b, respectively.

Analogously, Figure 5a,b illustrate the bending mode shapes of the same truncated AFGM beams, constructed based on the FEM results, which are given in Table 11.

It is worth noticing that the geometric non-uniformity and the variation of material parameters along the beam length significantly affect the mode shapes. The figures clearly show that the positions of the nodal points shift along the beam length, and this shift is more pronounced as the truncation factor ${\delta }_1$ becomes smaller.

4.5. AFGM Timoshenko Beams with Curvilinearly Double-Tapered Cross-Sections

To the best of authors' knowledge, benchmark results for curvilinearly double-tapered AFGM beams are not available in the existing literature. Hence, the results of the present analysis are validated against axial plumb and slender beams with Type 2 and Type 3 material distributions, as considered in [29], which feature variable thickness and constant width. Further, the analysis will be etended to include beams with double-tapered cross-sectional variations. The geometrical and material parameters are identical to those used in the previous examples of curvilinear taper beams with exponential and power-law gradations as presented in [29].

Figure 5. Mode shapes of truncated AFGM beams with parabolic taper height and constant width obtained from 3- D FEM analysis: (a) C-C boundary conditions; (b) C-F boundary conditions.

Figure 5. Mode shapes of truncated AFGM beams with parabolic taper height and constant width obtained from 3- D FEM analysis: (a) C-C boundary conditions; (b) C-F boundary conditions.

The first two dimensionless frequencies $\overline{\omega }=\omega L\sqrt{{\displaystyle \frac{{\rho }_{Al}}{E_{Al}}}}$ of axially graded Type 2 plump and slender C-C beams with variable thickness and constant width, as well beams with simultaneously variable thickness and width, depending on the exponential law parameter $a$ and gradient index ${\gamma }_e$ are collected in Table 12. Satisfactory agreement between the reference and present results for plump and slender AFGM beams with variable thickness is observed. The errors between the frequencies fall within the range observed for frequencies of similar beams with $a=1$ and ${\gamma }_e=1$ in Table 10. In particular, minor deviations are noted for the first frequency, while higher discrepancies are observed for the second frequency.

Upon reviewing Table 12, it becomes evident that introducing variations in width influences the frequencies of the beams. Specifically, plump double-tapered beams show slightly lower frequencies compared to those with only variable thickness, whereas slender double-tapered beams exhibit a slight increase. However, the trends in frequency changes with increasing gradient index ${\gamma }_e$ and parameter $a$ remain consistent for both beams with only variable thickness and double-tapered beams. The frequencies increase with higher ${\gamma }_e$ and $a$. This increase is more pronounced as ${\gamma }_e$ or $a$ become larger.

Analogously, axially graded Type 3 plump and slender C-C beams with variable thickness and constant width, and with simultaneously variable thickness and width are examined. Table 13 presents the first two dimensionless frequencies $\overline{\omega }=\omega L\sqrt{{\displaystyle \frac{{\rho }_{sus}}{E_{sus}}}}$ of these beams for different values of power law parameter $b$ and gradient index ${\gamma }_c$.

It is found in Table 13 that the agreement between the results obtained from the IGA analysis in [29] and the present study is satisfactory, with errors comparable to those observed for Type 3 plump and slender C-C beams with $b=1$ and ${\gamma }_c=1$ in Table 11. Additionally, Table 13 shows that introducing variable width along with variable height has similar effects on the frequencies of both beams with variable thickness only and curvilinearly double-tapered beams, similar to what was observed for Type 2 plump and slender C-C beams, as mentioned earlier.

On the other hand, the trends of the frequency changes with increasing $b$ and ${\gamma }_c$ are found to be opposite to those in Type 2 plump and slender C-C beams. In Type 3 plump and slender beams, the frequencies decrease with higher values of ${\gamma }_c$ and $b$. This decrease begins more rapidly and exhibits a faster decay as ${\gamma }_c$ or $b$ increase.

Table 14 presents the dimensionless natural frequencies of beams with variable thickness only and curvilinearly double-tapered beams under different boundary conditions. The calculation errors for the first and second frequencies are within 3%. Notably, the second bending frequency is not available from the IGA analysis in [29] in some cases. Analyzing the data in Table 14, it is evident that the frequencies of curvilinearly double-tapered beams behave similarly to those of beams with variable thickness only, depending on the boundary conditions. In particular, the highest frequencies correspond to C-C restraints, while the lowest frequencies are observed under C-F boundary conditions.

5. Conclusions

Free vibration analysis of axially functionally graded beams with curvilinear taper cross-sections and various material distributions is performed using the semi-analytical approach of the differential transform method. Unlike standard procedure that impose boundary conditions on transformed images of unknown functions, the present DTM approach integrates various boundary conditions directly into the original continuous problem formulation. The non-uniform beam shapes, defined as conic sections, are analytically described using rational Bézier curves in the Gröbner basis in a unified manner.

The equations of motion for free vibration are derived by principle of virtual work under the assumptions of the Timoshenko beam theory. It is assumed that both the material properties of the beam and its geometric features vary continuously along the axial direction.

Comparisons are presented to validate the accuracy and assess the efficiency of the proposed computational technique. The computed examples encompass both well-documented cases and previously unexamined scenarios, ranging from uniform homogeneous and inhomogeneous beams to curvilinearly double-tapered axially functionally graded beams. In the cases, where benchmark results are unavailable in the literature, the current DTM solutions are compared with our own FEM simulations. The three-dimensional finite element modeling was conducted using the ABAQUS code with a user-defined subroutine (UMAT) to incorporate material gradients into the 3-D element C3D20.

Based on these results, the DTM in its proposed formulation exhibits excellent performance in modeling the free vibration of beams with complex geometry and material distribution. This study provides new results from frequency analysis of curvilinearly double-tapered AFGM beams. It was found that introducing simultaneous variations in height and width of the beam cross-section influences the beam’s frequencies. While these changes may not be substantial, conducting accurate dynamic structural modeling remains crucial. Therefore, these findings can be useful in the design process of axially functionally graded beams and serve as benchmarks for researchers interested in exploring similar problems using alternative methods.