Free Vibration Analysis of Functionally Graded Carbon Nanotube-Reinforced Higher Order Refined Composite Beams Using Differential Quadrature Finite Element Method

Ihab Eddine Houalef, Ismail Bensaid, Ahmed Saimi^* and Abdelmadjid Cheikh

IS2M Laboratory, Faculty of Technology, Mechanical engineering Department, University, Abou Beckr Beklaid (UABT), Tlemcen, Algeria
E-mail: ahmedsaimi@hotmail.fr
*Corresponding Author

Received 01 March 2022; Accepted 28 November 2022; Publication 06 February 2023

Abstract

Present paper deals on the free vibration investigation of carbon nanotube-reinforced composite (CNTs) beams, based on refined third order shear deformation finite element beam theory. The particularity of this model is that, it can capture shear deformation effect without using of any shear correction factor by satisfying shear stress free at free edges. The carbon nanotubes are supposed to be immersed in a polymeric matrix with functionally graded pattern across the thickness direction of the beam, and their material properties are evaluated using the rule of mixture. The differential equations of motion and related boundary conditions are extracted using Lagrange’s principle and solved employing a robust numerical tool called, Differential Quadrature Finite Element Method (DQFEM) for the first time, with high convergence speed, fast calculus performance as well as a good numerical stability. The obtained results have been validated with those available in literature, in order to show the correctness of the present model. Afterwards, a deep parametric study is performed to examine the effects of various geometrical and material parameters on the vibration behavior of FG-CNTs beams.

Keywords: FG-CNTs beam, dynamic analysis, refined third order theory, differential quadrature finite elements method, enriched beam element.

1 Introduction

In modern engineering sectors, the use of fiber reinforced composite (FRC) to reinforce structural elements such as beams and plates has become a necessity, given the paramount importance that they provide, such as outstanding mechanical performance, sustainability, optical and low density. Among the new reinforcement methods are nano-composite materials, such as carbon nanotube-reinforced (CNTs) (Ijima 1991; Thostenson et al. 2001), which are considered as a new emerging nano composite materials because of the extreme benefit when using them as reinforcing elements, providing exceptional mechanical, thermal, electrical properties, and size dependency Kiang et al. (1998). To deal with several problems such as, delamination in laminated composite structures especially when using CNTs uniformly dispersed in the polymer matrix, the concept of functionally graded carbon nanotubes (FG-CTNs) has been proposed by Shen (2009) for designing the distributions of GPLs immersed in a polymer matrix by providing more performance. Based on this idea, various researches have been provided in the literature. Ke et al. (Ke et al. 2010, Ke et al. 2013) investigated the FGCNT volume fraction effects on dynamic stability of nonlinear vibration and of composite beams. Wang and Shen (Wang and Shen 2011) showed that the CNTRC plates with symmetrical form of CNTs provide lower natural frequencies compared to unsymmetrical or uniform distribution of CNTs, when the FG-CNTs plates are in a thermal environment. Wattanasakulpong and Ungbhakorn (2013) presented a deep study on the bending, buckling and vibration behaviors of carbon nanotube reinforced composite (CNTRC) beams, considering various higher-order shear deformation beam theories. Shen and Xiang (2013) investigated the thermal postbuckling response of FG-CNTRC beams considering both edges simply supported, employing perturbation method. The element-free kp-Ritz method has been used by Lei et al. (2013), to investigate the free vibration of FG-CNT reinforced composite rectangular plates in a thermal environment. Lin and Xiang (2014) used both first and third order beam theories in conjunction with p-Ritz method, for investigating free vibration response of aligned and graded CNTs composite beams. Yang et al. 2015 researched the dynamic buckling behaviors of FG nanocomposite beams strengthened by CNT as a core and integrated with two surface bonded piezoelectric layers. Wu et al. 2015 examined the free vibration and buckling performance of sandwich beams strengthened with FGCNTRCs face sheets employing the Timoshenko beam model. Tagrara et al. 2015 provided an original trigonometric beam theory to study the bending, buckling and vibration responses of functionally graded carbon nanotube-reinforced composite beams. The shear buckling behaviour of FG-CNTRC composite plates has been explored by Kiani Y (2016), taking into account various types of boundary conditions employing the Chebychev–Ritz technique. Jam and Kiani (2015) studied the behavior of FG-CNTRC reinforced beam subjected to a low velocity impact due to a single mass in thermal environment. Ebrahimi and Karimiasl (2018) presented an analytical approach to explore the surface and flexoelectric effects on the buckling characteristics of an embedded piezoelectric sandwich nanobeam. Kiani Y (2018) examined the thermal post-buckling response of a sandwich beam made with CNTRC face sheets and resting on a Pasternak foundation. They showed that the properties of CNTs are important factors on thermal buckling and postbuckling characteristics of the sandwich beam with FG-CNTRC face sheets. VO-DUY et al. (2019) employed a linear two-node element with six degrees of freedom combined with the first-order shear deformation theory, to study the free vibration of laminated functionally graded carbon nanotube reinforced composite beams. Babaei et al. (2021) contributed on the vibrational behavior investigation of thermally pre-/post buckled FG-CNTRC beams resting on a nonlinear elastic foundation using a two-step perturbation procedure.

Numerical and semi-numerical methods, such as the popular finite element method (FEM) and differential quadrature have found wide usage in recent decades, due to their efficiency, flexibility and adaptability in solving differential equations related to structural elements with complicated effects, related to geometries, shapes, interactions, boundary conditions, etc. with precise results. Since then, several studies have been provided lately based on these numerical methods. Vo and Thai (2012) and Vo et al. (2014) developed a two-noded C${}^1$ beam element with five degree-of-freedom per node, to investigate the static deflection, vibration and buckling behaviour of composite and functionally graded (FG) beams with refined shear deformation theory. In another study, a quasi-3D theory was applied by Vo et al. (2017) for studying the free vibration of axially loaded composite beams, based on two-noded C${}^1$ with six degrees of freedom. Yarasca et al. (2016) presented a Hermitian–Lagrangian finite element formulation of a generalized quasi-3D hybrid type HSDT considering 7DOF, for static analysis of functionally graded single and sandwich beams. Forced Vibration of was studied by Barati and Shahverdi (2020) employing the refined shear deformation beam theory in conjunction with two-noded C${}^1$ finite beam element, which contains ten degrees of freedom, for inspecting forced vibration of a nanocomposite beam reinforced with different distributions of graphene platelets (GPLs) in thermal environments. Karamanli and Vo (2021) presented Finite elements solutions based on normal and shear deformation beam theory, to analyze the bending, buckling and free vibration of CNTRC/GPLRC beams. Azimi et al. (2018) analyzed the vibration of rotating, functionally graded Timoshenko nano-beams under an in-plane nonlinear thermal loading, via the differential quadrature method. The GDQM was used by Lal and Dangi (2019), to explore the effect of linear and nonlinear thermal environment on vibration characteristic of temperature dependent BD-FG non-uniform nanobeam on the basis of Timoshenko beam theory together with Eringen’s nonlocal elasticity theory. Lei et al. (2019) researched the postbuckling response of bi-directional FG beams with porosities based on a novel third-order shear deformation theory. Even and uneven distributions of porosities were considered in their investigation. Abdollahi and Yas (2020) studied free vibration response of boron nitride nanotubes (BNNTs) reinforced Timoshenko beams on an elastic foundation. The orientation of nanotubes in the matrix was considered to be unidirectional or randomly oriented. Natural frequencies were gained by introducing (GDQM). Bensaid and Saimi (2022) explored the dynamic behavior of functionally graded porous beams resting on viscoelastic foundation, employing the generalized differential quadrature method. Eltaher et al. (2020) performed a study on the static stability response and mode-shapes of laminated composite beams and undergoing varying axial inplane loads with differential quadrature method (DQM). The nonlinear thermal buckling and postbuckling of bidirectional functionally graded non-uniform microbeams were inspected by Attia and Mohamed (2020) employing Reddy beam model and generalized differential quadrature method (GDQM). It is well known that the coupling between numerical methods has an extreme advantage, and it has attracted the attention of several researchers thanks to the great advantages it provides such as, improved numerical stability, high computation efficiency, fast convergence, complex geometries, etc. Besides, works on the combination between the finite element method and the generalized differential quadrature method is somewhat limited in the literature. One can cite among the first works provided by Xing and Liu, (2009), which presented a Differential Quadrature Finite Element Method (DQFEM) to facilitate the complexity of imposing boundary conditions in DQM with high accuracy and results convergence in application to free vibrations of thin plate with curvilinear domain. Liu C et al. (2016) proposed a differential quadrature hierarchical finite element method (DQHFEM) that introduce interpolation basis on the boundary of hierarchical finite element method elements, to study the vibration and bending of Mindlin plates with curvilinear domains. Yan et al. (2021) developed a novel unified quasi-3D solution based on combination between Carrera Unified Formulation (CUF) and differential quadrature finite element method (DQFEM), to explore Free vibration analysis of variable stiffness composite laminated beams and plates. Recently, Saimi et al. (2021) extended the differential quadrature finite element method (DQFEM) and hierarchical finite element method (DQHFEM) to examine the dynamic behavior of on-board shaft, for the first time. As we can see previously, there was no reported work done on free vibration of refined higher order shear deformation beams reinforced with a functionally graded carbon nanotube using (FGCNTs). This paper aims to investigate for the first time within the framework of refined higher order shear deformation beam theory (RSHBT) strengthened with functionally graded carbon nanotubes. The material properties of carbon nanotubes are assumed to vary in the thickness direction in an FG pattern. Differential quadrature coupled with finite element method (DQFEM) as a power numerical tool is employed to solve the differential governing equations of composite beams for the first time. After that, a parametric investigation is conducted to explore the impacts of length thickness ratios, carbon nanotube volume fraction and their distribution forms, combined physical and geometrical parameters in three-dimension, various boundary conditions on the free vibration behaviors of FG-CNTRC composite beams.

2 FG-CNRC Beams

In this study, we consider a straight composite beam made with a mixture of isotropic polymer matrix and SWCNTs (10,10) armchair dispersed in it. Figure 1(a) shows a CNTRC beam, having length (L) and thickness (h). In the present investigation, four different patterns of reinforcement through the cross sections are considered as shown in Figure 1(b).

Figure 1(a): A geometry of FG-CNTs composite beam.

Figure 1(b): Different patterns of CNTs reinforcement through a beam section.

The material properties of CNTRC beams are evaluated employing the extended rule of mixture which provides the effective Young’s modulus and shear modulus of CNTRC beams as (Shen 2009, Wattanasakulpong and Ungbhakorn 2013).

$$E_{11}={\eta }_1{V}_{cnt}{E}_{11}^{cnt}+V_p{E}^p$$	(1a)
$$\frac{{\eta }_2}{E_{22}}=\frac{V_{cnt}}{E_{22}^{cnt}}+\frac{V_p}{E^p}$$	(1b)
$$\frac{{\eta }_3}{G_{12}}=\frac{V_{cnt}}{G_{12}^{cnt}}+\frac{V_p}{G^p}$$	(1c)

in which, $E_{11}^{\mathit{\text{cnt}}}$; $E_{22}^{\mathit{\text{cnt}}}$ and $G_{12}^{\mathit{\text{cnt}}}$ represent the Young’s modulus and shear modulus of SWCNT, correspondingly and E${}^p$ and G${}^p$ are the corresponding material properties of the polymer matrix. Also, V${}_{c\mathit{}n\mathit{}t}$ and V${}_p$ are the volume fractions for carbon nanotube and the polymer matrix, respectively, with the relation of $V_{cnt}+V_p=1$. To consider the size-dependency material properties of SWCNT, the CNT efficiency parameters, ${\eta }_i$ (i $=$ 1, 2, 3), are selected. Theses parameters have been obtained by matching the elastic moduli of CNTRCs assessed by the MD simulation by the numerical results determined by the rule of mixture (Han and Elliott 2007). Utilizing the same rule, Poisson’s ratio (v) and mass density ($\rho$) of the CNTRC beams are expressed as

$$\nu =V_{cnt}{\nu }^{cnt}+V_p{\nu }^p;\rho =V_{cnt}{\rho }^{cnt}+V_p{\rho }^p$$

(2)

where ${\nu }^{cnt}$, ${\nu }^p$ and ${\rho }^{cnt}$, ${\rho }^p$ are the Poisson’s ratios and densities of the CNT and polymer matrix respectively. Explicit mathematical term of CNTs volume fraction in every case of distribution is given in Table 1.

Table 1 Volume fraction of CNTs as a function of thickness coordinate for various cases of CNTs distribution

CNTs distribution	V${}_{c\mathit{}n}$
UD CNTRC	$V_{cn}^{*}$
FG-V CNTRC	$V_{cn}^{*}\left(1+2{\displaystyle \frac{z}{h}}\right)$
FG-O CNTRC	$2V_{cn}^{*}\left(1-2{\displaystyle \frac{\left|z\right|}{h}}\right)$
FG-X CNTRC	$4V_{cn}^{*}{\displaystyle \frac{\left|z\right|}{h}}$

where $V_{cnt}^{*}$ represents the volume fraction of CNTs, which can be evaluated from thefollowing equation

$$V_{cnt}^{*}=\frac{W_{cnt}}{W_{cnt}+({\rho }^{cnt}/{\rho }^m)(1-W_{cnt})}$$

(3)

In which W${}_{c\mathit{}n\mathit{}t}$ is the mass fraction of CNTs related to the composite beam. From the above table, one can observe that the V-, O-and X-Beams are sorts of functionally graded shapes in which their material components are changed continuously across their thicknesses; while, the UD-Beam has uniformly distributed CNT reinforcement. The parameter ($\eta$${}_i$) cited previously has the following numerical values ${\eta }_1$ $=$ 1.2833 and ${\eta }_2$ $=$ ${\eta }_3$ $=$ 1.0566 for the case of $V_{cnt}^{*}$ ;$=$ 0.12; ${\eta }_1$ $=$ 34.14 and ${\eta }_2$ $=$ ${\eta }_3$ $=$ 1.7101 for the case of $V_{cnt}^{*}$ ;$=$ 0.17; ${\eta }_1$ $=$ 1.3238 and ${\eta }_2$ $=$ ${\eta }_3$ $=$ 1.7380 for the case of $V_{cnt}^{*}$ $=$ 0.28 (Yas and Samadi 2012).

3 Basic Mathematical Modeling

3.1 Kinematics and Equations of Motion

The general equations of motion related to this investigation for investigating the dynamic behavior of FG-CNTs composite beam are extracted from refined shear deformation beam theory (RSDBT) without the need to use a shear correction factor (SCF). The displacement field of any arbitrary point in the mid plan of the beams alongside the x- and z-axes can be given as follows (Thai and Vo 2012; Vo et al. 2014; Bekhadda et al. 2019):

$$u(x,z,t)=u_0(x,t)-z\frac{d{w}_b}{dx}-f(z)\frac{d{w}_s}{dx}$$	(4a)
$$v(x,z,t)=0$$	(4b)
$$w(x,z,t)=w_b(x,t)+w_s(x,t)$$	(4c)

By which, $u_0(x,t)$, $w_b(x,t)$ and $w_s(x,t)$ are successively, the in-plane displacement in x-ways, bending and shear components of the transverse displacement of points on the neutral axis of the beam; and $f(z)$ is a shape function determining the distribution of the transverse shear strain and shear stress over the depth of the composite beam.

The expression of normal and shear strain components associated with the displacement field in Equation (2) are given by

$${\epsilon }_{xx}=\frac{d{u}_0}{dx}-z\frac{d^2{w}_b}{d{x}^2}-f(z)\frac{d^2{w}_s}{d{x}^2}$$
$${\gamma }_{xz}=g(z)\frac{d{w}_s}{dx}$$	(5)
$$\text{Where}f(z)=z-\frac{h}{dx}\sin \frac{\pi z}{h},g(z)=1-\frac{df(z)}{dz}.$$

By assuming that the material of CNTs beam follows Hooke’s law, the stresses in the composite beam become

$$Q_{11}(z)=\frac{E_{11}(z)}{1-{\nu }^2}$$	(6a)
$$Q_{11}(z)=G_{12}(z)$$	(6b)

3.2 Equations of Motion

In order to derive the general equations of motion, Lagrange’s principle is used here in (Saimi, Bensaid et al. 2021) as follows:

$$\frac{d}{dt}\left(\frac{\partial \mathrm{\Pi }}{\partial \dot{q_i}}\right)+\frac{\partial \mathrm{\Pi }}{\partial q_i}=0$$

(7)

where q${}_i$represents the unknown coefficients (U${}_i$, W${}_{b\mathit{}i}$ and W${}_{s\mathit{}i}$), and ($\mathrm{\Pi }$) is total energy functional, the over-dot denotes the partial derivative with respect to time.

The variation of the strain energy of the beam is of the form

$U_e^B$	$={\displaystyle \frac{1}{2}}{\displaystyle {\int }_{V_e}}{\left\{{\displaystyle \genfrac{}{}{0.0pt}{}{{\sigma }_x}{{\tau }_{xz}}}\right\}}^T\left\{{\displaystyle \genfrac{}{}{0.0pt}{}{{\epsilon }_x}{{\gamma }_{xz}}}\right\}\mathit{d}V$
$U_e^B$	$={\displaystyle \frac{1}{2}}{\displaystyle {\int }_{V_e}}{\left\{{\displaystyle \genfrac{}{}{0.0pt}{}{{\epsilon }_x}{\tau {\gamma }_{xz}}}\right\}}^T\left[{\displaystyle \genfrac{}{}{0.0pt}{}{E}{0}}{\displaystyle \genfrac{}{}{0.0pt}{}{\mathrm{\hspace{0.25em}\hspace{0.25em}\hspace{0.25em}\hspace{0.25em}\hspace{0.25em}\hspace{0.25em}0}}{G}}\right]\left\{{\displaystyle \genfrac{}{}{0.0pt}{}{{\epsilon }_x}{{\gamma }_{xz}}}\right\}\mathit{d}V$
$U_e^B$	$={\displaystyle \frac{1}{2}}b{\displaystyle {\int }_0^l}\left\{{\left[\begin{array}{c}\hfill {\displaystyle \frac{d{u}_0}{dx}}\hfill \\ \hfill {\displaystyle \frac{d^2{w}_b}{d{x}^2}}\hfill \\ \hfill {\displaystyle \frac{d^2{w}_s}{d{x}^2}}\hfill \end{array}\right]}^T\left[\begin{array}{ccc}\hfill I_1\hfill & \hfill I_2\hfill & \hfill I_3\hfill \\ \hfill I_2\hfill & \hfill I_5\hfill & \hfill I_4\hfill \\ \hfill I_3\hfill & \hfill I_4\hfill & \hfill I_6\hfill \end{array}\right]\left[\begin{array}{c}\hfill {\displaystyle \frac{d{u}_0}{dx}}\hfill \\ \hfill {\displaystyle \frac{d^2{w}_b}{d{x}^2}}\hfill \\ \hfill {\displaystyle \frac{d^2{w}_s}{d{x}^2}}\hfill \end{array}\right]+{\displaystyle \frac{d{w}_s}{dx}}{I}_7{\displaystyle \frac{d{w}_s}{dx}}\right\}\mathit{d}x$	(8)
$U_e{}^B$	$={\displaystyle \frac{1}{2}}b{\displaystyle {\int }_0^l}[(I_1{\displaystyle \frac{d^2{u}_0}{d{x}^2}}+2I_2{\displaystyle \frac{d{u}_0}{dx}}{\displaystyle \frac{d^2{w}_b}{d{x}^2}}+2I_3{\displaystyle \frac{d{u}_0}{dx}}{\displaystyle \frac{d^2{w}_s}{d{x}^2}}+2I_4{\displaystyle \frac{d^2{w}_b}{d{x}^2}}{\displaystyle \frac{d^2{w}_s}{d{x}^2}}$
	$\mathrm{\hspace{1em}}+I_5{\displaystyle \frac{d^2{w}_b}{d{x}^2}}{\displaystyle \frac{d^2{w}_b}{d{x}^2}}+I_6{\displaystyle \frac{d^2{w}_s}{d{x}^2}}{\displaystyle \frac{d^2{w}_s}{d{x}^2}}+I_7{\displaystyle \frac{d{w}_s}{dx}}{\displaystyle \frac{d{w}_s}{dx}})]dx$	(9)

where $V_e$ is the volume of the element; $I_1,I_2,I_3,I_4,I_5,I_6,I_7$ are the beam rigidities, defined as

$$(I_1,I_2,I_3,I_4,I_5,I_6)=b{\int }_{-h/2}^{h/2}{Q}_{11}(1,z,f,zf,z^2,f^2)$$	(10)
$$I_7=b{\int }_{-h/2}^{h/2}{Q}_{55}(g^2)$$	(11)

The variation of kinetic energy of the beam can be expressed as

$T_e$	$={\displaystyle \frac{1}{2}}{\displaystyle {\int }_{V_e}}{\rho }_f(\dot{u^2}+\dot{w^2})\mathit{d}V$
$T_e$	$={\displaystyle \frac{1}{2}}{\displaystyle {\int }_{V_e}}{\rho }_f\left[{\left\{{\dot{u}}_0-z{\dot{w}}_b-f(z){\displaystyle \frac{d{\dot{w}}_s}{dx}}\right\}}^2+\left(\dot{w}_b{}^2+\dot{w}_s{}^2\right)\right]\mathit{d}V$
$T_e$	$={\displaystyle \frac{1}{2}}b{\displaystyle {\int }_0^l}\left\{{\left[\begin{array}{c}\hfill {\dot{u}}_0\hfill \\ \hfill {\dot{w}}_b\hfill \\ \hfill {\dot{w}}_s\hfill \end{array}\right]}^T\left[\begin{array}{ccc}\hfill J_1\hfill & \hfill J_2\hfill & \hfill J_3\hfill \\ \hfill J_2\hfill & \hfill J_5\hfill & \hfill J_4\hfill \\ \hfill J_3\hfill & \hfill J_4\hfill & \hfill J_6\hfill \end{array}\right]\left[\begin{array}{c}\hfill {\dot{u}}_0\hfill \\ \hfill {\dot{w}}_b\hfill \\ \hfill {\dot{w}}_s\hfill \end{array}\right]+{\left[\begin{array}{c}\hfill {\dot{w}}_b\hfill \\ \hfill {\dot{w}}_s\hfill \end{array}\right]}^T\left[\begin{array}{cc}\hfill J_1\hfill & \hfill J_1\hfill \\ \hfill J_1\hfill & \hfill J_1\hfill \end{array}\right]\left[\begin{array}{c}\hfill {\dot{w}}_b\hfill \\ \hfill {\dot{w}}_s\hfill \end{array}\right]\right\}\mathit{d}x$	(12)
$T_e$	$={\displaystyle \frac{1}{2}}b{\displaystyle {\int }_0^l}[J_1({\dot{u^2}}_0+\dot{w}_b{}^2+\dot{w}_s{}^2+2{\dot{w}}_b{\dot{w}}_s)-2J_2{\dot{u}}_0{\displaystyle \frac{d{\dot{w}}_b}{dx}}-2J_3{\dot{u}}_0{\displaystyle \frac{d{\dot{w}}_s}{dx}}$
	$\mathrm{\hspace{1em}}+2J_4{\displaystyle \frac{d{\dot{w}}_b}{dx}}{\displaystyle \frac{d{\dot{w}}_s}{dx}}+J_5{\left({\displaystyle \frac{d{\dot{w}}_b}{dx}}\right)}^2+J_6{\left({\displaystyle \frac{d{\dot{w}}_s}{dx}}\right)}^2]dx$

where the over dot denotes the derivative with respect the time variable t, and the mass moment $J_1,J_2,J_3,J_4,J_5,J_6$s are defined as

$$(J_1,J_2,J_3,J_4,J_5,J_6)=b{\int }_{-h/2}^{h/2}\rho (1,z,f,zf,z^2,f^2)$$

(14)

3.3 Solution Procedure Based DQFEM

3.3.1 The reformulated differential quadrature rule

Differential quadrature rules approximate the derivatives of a function using a weighted linear sum of field variables along a line passing through the point. For polynomial basis functions DQ, a set of Lagrange polynomials are employed as the test functions.

Thus, for a field variable $f(x)$ its derivative of order n in a discrete point $x_i$ can be expressed as:

$$\frac{{\partial }^{n}f(x;t)}{\partial x^n}=\sum _{j=1}^N{A}_{ij}^{(n)}f(x_j,t)\mathit{\hspace{1em}}(i=1,2,3,\mathrm{\dots },N)$$

(15)

Where $A_{ij}^{(n)}$ is the weighting coefficient related to the derivative of order n, and the weighting coefficient is obtained as follows if $n=1$, so

$$A_{ij}^{(1)}=\frac{M(x_i)}{(x_i-x_j)M(x_j)}\mathit{\hspace{1em}}i\ne j,j=1,2,\mathrm{\dots },N$$
$$A_{ii}^{(1)}=-\sum _{j=1,j\ne i}^n{A}_{ij}^{(1)}\mathit{\hspace{1em}}i=1,2,\mathrm{\dots },N$$	(16)

where

$M(x_i)$	$={\displaystyle \prod _{k=1,k\ne i}^N}(x_i-x_k)$
$M(x_j)$	$={\displaystyle \prod _{k=1,k\ne i}^N}(x_j-x_k)$	(17)

If n $>$ 1, secondary and higher order derivatives, the weighting coefficients are determined using the following simple recurrence relationship:

$$A_{ij}^{(n)}=n\left(A_{ij}^{(1)}*A_{ii}^{(n-1)}-\frac{A_{ij}^{(n-1)}}{(x_i-x_j)}\right)\mathit{\hspace{1em}}i\ne j,i,j=1,2,\mathrm{\dots },N,n>1$$
$$A_{ii}^{(n)}=-\sum _{j=1,j\ne i}^N{A}_{ij}^{(n)}\mathit{\hspace{1em}}i=1,2,\mathrm{\dots },N$$	(18)

3.3.2 Gauss-Lobatto quadrature rule

The theory of Gauss-Lobatto quadrature rules can be found in the mathematical literature; The Gauss Lobatto quadrature rule with a degree of accuracy (2n$-$3) for the function f(x) defined in [$-$1; 1] is:

$${\int }_{-1}^{1}f(x)\mathit{d}x=\sum _{j=1}^N{C}_{j}f(x_j)$$

(19)

With the weighting coefficient $C_j$ of the Gauss-Lobatto, integration is given by:

$$C_1=C_N=\frac{2}{N(N-1)},C_j=\frac{2}{N(N-1){[P_{N-1}(x_j)]}^2}\mathit{\hspace{1em}}(j\ne 1,N)$$

(20)

$x_j$ is the $(j-1)$ zero of the first order derivative of $P_{N-1}(x)$. To solve the roots of the Legendre polynomials, we will use the recursivity formula as Equations (21) and (22), it is easy to obtain thousands of roots.

$$P_{N+1}(x)=\frac{2N+1}{N+1}x{P}_N(x)-\frac{N}{N+1}{P}_{N-1}(x)$$

(21)

With $P_0(x)=1$, $P_1(x)=x$. The n${}^{\mathrm{th}}$-order derivation of the Legendre polynomials can be determined by the following formula:

$$P_{N+1}^{(n)}(x)=x{P}_N^{(n)}(x)+(N+n)P_N^{(n)}(x)$$

(22)

In order to obtain a denser population near the boundaries, sampling points are selected according to the grid distribution of Gauss–Lobatto nodes.

$$x_j=-\cos \left(\frac{j-1}{N-1}\pi \right)$$

(23)

Gauss-Lobatto nodes are solved with the Newton-Raphson iteration method.

$$x^{iT+1}=x^{iT}-F^{\prime }{(x^{iT})}^{-1}F(x^{iT}),iT=0,1,\mathrm{\dots }$$

(24)

in which

$$x={[x_2,x_3,\mathrm{\dots },x_{N-1}]}^T$$	(25)
$$F(x)={[f(x_2),f(x_3),\mathrm{\dots },f(x_{N-1})]}^T$$	(26)
$$F^{\prime }(x)={\left[\frac{\partial f(x_j)}{\partial x_i}\right]}_{(N-2)\times (N-2)}$$	(27)
$$f(x_j)=\sum _{k-1,k\ne j}^N\frac{1}{x_j-x_k}\mathit{\hspace{1em}}j=2,3,\mathrm{\dots },N-1$$	(28)
$$\frac{\partial f(x_j)}{\partial x_i}=\{\begin{array}{cc}-\sum _{k=1,k\ne j}^N\frac{1}{{(x_j-x_k)}^2},\hfill & (i=j)\hfill \\ \frac{1}{{(x_j-x_k)}^2},\hfill & (i\ne j)\hfill \end{array}$$	(29)

Where k is the value of x at i${}^{{\mathrm{t}}^{\mathrm{th}}}$ iteration step. This method is less sensitive to the initial value. The values given by Equation (29) are used as initial values.

3.3.3 The differential quadrature finite element method

The finite element method (FEM) is considered as an efficient tool for the numerical solution, which is used in wide range of engineering problems. However, this method sometimes lacks in both convergence and calculation speed and numerical instability; currently it is necessary to couple this method with other numerical or semi-numerical methods, in order to make it more robust, stable and faster in the calculus. The differential quadrature finite element method was developed by (Xing and Liu, 2009), whose differential quadrature rules and Gauss-Lobatto quadrature are used to discretize the system energies.

Assuming that the deflection function is

$u(x)$	$={\displaystyle \sum _{i=1}^N}{L}_i(x)u_i$
$w(x)$	$={\displaystyle \sum _{i=1}^N}{L}_i(x)w_i$	(30)

With $L_i$ is the Lagrange polynomial, and $u_i=u(x_i)$, $w_i=w(x_i)$ are the displacements of the Gauss Lobatto quadrature points or the DQ nodal displacements of the beam finite element.

Using DQ rules and Gauss - Lobatto quadrature the expressions of kinetic energy and strain energy (8–12) can be written as follows:

$T_e$	$={\displaystyle \frac{1}{2}}b[J_1([{\overline{Q}}^{-T}\overline{C}{\overline{Q}}^{-1}]\dot{u}_0{}^2+[{\overline{Q}}^{-T}\overline{C}{\overline{Q}}^{-1}]\dot{w}_b{}^2+[{\overline{Q}}^{-T}\overline{C}{\overline{Q}}^{-1}]\dot{w}_s{}^2$
	$\mathrm{\hspace{1em}}+2[{\overline{Q}}^{-T}\overline{C}{\overline{Q}}^{-1}]{\dot{w}}_b{\dot{w}}_s)-2J{}_2[{\overline{Q}}^{-T}\overline{C}{\overline{A}}^{(1)}{\overline{Q}}^{-1}]\dot{u}{}_0\dot{w}{}_b$
	$\mathrm{\hspace{1em}}-2J_3[{\overline{Q}}^{-T}\overline{C}{\overline{A}}^{(1)}{\overline{Q}}^{-1}]{\dot{u}}_0{\dot{w}}_s+2J_4[{\overline{Q}}^{-T}{\overline{A}}^{(1)}\overline{C}{\overline{A}}^{(1)}{\overline{Q}}^{-1}]{\dot{w}}_s{\dot{w}}_b$
	$\mathrm{\hspace{1em}}+J_5[{\overline{Q}}^{-T}{\overline{A}}^{(1)}\overline{C}{\overline{A}}^{(1)}{\overline{Q}}^{-1}]\dot{w}_b{}^2+J_6[{\overline{Q}}^{-T}{\overline{A}}^{(1)}\overline{C}{\overline{A}}^{(1)}{\overline{Q}}^{-1}]\dot{w}_s{}^2]$	(31)
$U_e{}^B$	$={\displaystyle \frac{1}{2}}b{I}_1[u_0{}^T{\overline{Q}}^{-T}\overline{A}^{(1)}{}^T\overline{C}{\overline{A}}^{(1)}{\overline{Q}}^{-1}{u}_0]$
	$\mathrm{\hspace{1em}}+2I_2[u_0{}^T{\overline{Q}}^{-T}\overline{A}^{(1)}{}^T\overline{C}{\overline{A}}^{(2)}{\overline{Q}}^{-1}{w}_b]$
	$\mathrm{\hspace{1em}}+2I_3[u_0{}^T{\overline{Q}}^{-T}\overline{A}^{(1)}{}^T\overline{C}{\overline{A}}^{(2)}{\overline{Q}}^{-1}{w}_s]$
	$\mathrm{\hspace{1em}}+2I_4[w_b{}^T{\overline{Q}}^{-T}\overline{A}^{(2)}{}^T\overline{C}{\overline{A}}^{(2)}]$
	$\mathrm{\hspace{1em}}+2I_4[w_b{}^T{\overline{Q}}^{-T}\overline{A}^{(2)}{}^T\overline{C}{\overline{A}}^{(2)}{\overline{Q}}^{-1}{w}_s]$
	$\mathrm{\hspace{1em}}+I_5[w_b{}^T{\overline{Q}}^{-T}\overline{A}^{(2)}{}^T\overline{C}{\overline{A}}^{(2)}{\overline{Q}}^{-1}{w}_b]$
	$\mathrm{\hspace{1em}}+I_6[w_s{}^T{\overline{Q}}^{-T}\overline{A}^{(2)}{}^T\overline{C}{\overline{A}}^{(2)}{\overline{Q}}^{-1}{w}_s]$
	$\mathrm{\hspace{1em}}+I_7[w_s{}^T{\overline{Q}}^{-T}\overline{A}^{(1)}{}^T\overline{C}{\overline{A}}^{(1)}{\overline{Q}}^{-1}{w}_s]$	(32)

With $A^{(1)}$and $A^{(2)}$ indicates the matrices of the weighting coefficients of the DQ rules for the first and second order derivatives respectively calculated with Equations (16–18), with respect to the Gauss Lobatto nodes, and

$$C=\mathit{diag}[C_1{C}_2\mathrm{\dots }{C}_N]$$

(33)

Where $C_j$ are the weighting coefficients of the Gauss-Lobatto integration.

${\overline{u}}^T$	$=[u_1{u}_2\mathrm{\dots }{u}_N]$
${\overline{w}}^T$	$=[w_1{w}_2\mathrm{\dots }{w}_N]$	(34)

In order to construct an element that satisfies the requirements of continuity between elements, the element displacement vectors must be:

$u^T$	$=[u_1{u}_3^{\prime }{u}_3\mathrm{\dots }{u}_{N-2}{u}_N{u}_N^{\prime }]$
$w^T$	$=[w_1{w}_3^{\prime }{w}_3\mathrm{\dots }{w}_{N-2}{w}_N{w}_N^{\prime }]$	(35)

The relation between u and w is defined using the DQ rule:

$$u=Q\overline{u},w=Q\overline{w}$$

(36)

Where

$$Q=\left[\begin{array}{cccccc}\hfill 1\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill \mathrm{\cdots }\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill A_{1,1}^{(1)}\hfill & \hfill A_{1,2}^{(1)}\hfill & \hfill A_{1,3}^{(1)}\hfill & \hfill \mathrm{\cdots }\hfill & \hfill A_{1,N-1}^{(1)}\hfill & \hfill A_{1,N}^{(1)}\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill 1\hfill & \hfill \mathrm{\cdots }\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill \mathrm{⋮}\hfill & \hfill \mathrm{⋮}\hfill & \hfill \mathrm{⋮}\hfill & \hfill \mathrm{\ddots }\hfill & \hfill \mathrm{⋮}\hfill & \hfill \mathrm{⋮}\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill \mathrm{\cdots }\hfill & \hfill 0\hfill & \hfill 1\hfill \\ \hfill A_{N,1}^{(1)}\hfill & \hfill A_{N,2}^{(1)}\hfill & \hfill A_{N,3}^{(1)}\hfill & \hfill \mathrm{\cdots }\hfill & \hfill A_{N,N-1}^{(1)}\hfill & \hfill A_{1,N}^{(1)}\hfill \end{array}\right]$$

(37)

All forms of node distribution for differentiation and quadrature are [$-$1; 1]. Therefore, in order to apply them in practice, the following modifications must be made to the differential and quadrature matrices,

$$\overline{C}=\frac{l_e}{2}C,{\overline{A}}^{(1)}=\frac{2}{l_e}{A}^{(1)},{\overline{A}}^{(2)}=\frac{4}{l_e}{A}^{(2)}$$

(38)

Where $l_e$ is the length of the beam element.

The matrices for the entire system are obtained according to the FEM rules for assembling elementary matrices,

$$[M]\left\{\begin{array}{c}\hfill {\ddot{u}}_0(t)\hfill \\ \hfill {\ddot{w}}_b(t)\hfill \\ \hfill {\ddot{w}}_s(t)\hfill \end{array}\right\}+[K]\left\{\begin{array}{c}\hfill u_0(t)\hfill \\ \hfill w_b(t)\hfill \\ \hfill w_s(t)\hfill \end{array}\right\}=[0]$$

(39)

The elementary mass matrix obtained with (DQFEM)

$[M]=\left[\begin{array}{ccc}{J}_1\left[{\overline{Q}}^{-T}\overline{C}{\overline{Q}}^{-1}\right]\mathit{\hspace{1em}}\hfill & -J_2\left[{\overline{Q}}^{-T}\overline{C}{\overline{A}}^{\left(1\right)}{\overline{Q}}^{-1}\right]\mathit{\hspace{1em}}\hfill & -J_3\left[{\overline{Q}}^{-T}\overline{C}{\overline{A}}^{\left(1\right)}{\overline{Q}}^{-1}\right]\hfill \\ -J_2\left[{\overline{Q}}^{-T}\overline{C}{\overline{A}}^{\left(1\right)}{\overline{Q}}^{-1}\right]\mathit{\hspace{1em}}\hfill & J_1\left[{\overline{Q}}^{-T}\overline{C}{\overline{Q}}^{-1}\right]\mathit{\hspace{1em}}\hfill & J_1\left[{\overline{Q}}^{-T}\overline{C}{\overline{Q}}^{-1}\right]\hfill \\ & +J_5[{\overline{Q}}^{-T}\overline{A}^{\left(1\right)}{}^T\mathit{\hspace{1em}}\hfill & +J_4[{\overline{Q}}^{-T}\overline{A}^{\left(1\right)}{}^T\hfill \\ & \overline{C}{\overline{A}}^{\left(1\right)}{\overline{Q}}^{-1}]\mathit{\hspace{1em}}\hfill & \overline{C}{\overline{A}}^{\left(1\right)}{\overline{Q}}^{-1}]\hfill \\ -J_3\left[{\overline{Q}}^{-T}\overline{C}{\overline{A}}^{\left(1\right)}{\overline{Q}}^{-1}\right]\mathit{\hspace{1em}}\hfill & J_1\left[{\overline{Q}}^{-T}\overline{C}{\overline{Q}}^{-1}\right]\mathit{\hspace{1em}}\hfill & J_1\left[{\overline{Q}}^{-T}\overline{C}{\overline{Q}}^{-1}\right]\hfill \\ & +J_4[{\overline{Q}}^{-T}\overline{A}^{\left(1\right)}{}^T\mathit{\hspace{1em}}\hfill & +J_6[{\overline{Q}}^{-T}\overline{A}^{\left(1\right)}{}^T\hfill \\ & \overline{C}{\overline{A}}^{\left(1\right)}{\overline{Q}}^{-1}]\mathit{\hspace{1em}}\hfill & \overline{C}{\overline{A}}^{\left(1\right)}{\overline{Q}}^{-1}]\hfill \end{array}\right]$

(40)

The stiffness matrices obtained with (DQFEM)

$[K]=\left[\begin{array}{ccc}{I}_1\left[{\overline{Q}}^{-T}\overline{A}^{\left(1\right)}{}^T\overline{C}{\overline{A}}^{\left(1\right)}{\overline{Q}}^{-1}\right]\mathit{\hspace{1em}}\hfill & I_2\left[{\overline{Q}}^{-T}\overline{A}^{\left(1\right)}{}^T\overline{C}{\overline{A}}^{\left(2\right)}{\overline{Q}}^{-1}\right]\mathit{\hspace{1em}}\hfill & I_3\left[{\overline{Q}}^{-T}\overline{A}^{\left(1\right)}{}^T\overline{C}{\overline{A}}^{\left(2\right)}{\overline{Q}}^{-1}\right]\hfill \\ I_2\left[{\overline{Q}}^{-T}\overline{A}^{\left(2\right)}{}^T\overline{C}{\overline{A}}^{\left(1\right)}{\overline{Q}}^{-1}\right]\mathit{\hspace{1em}}\hfill & I_5\left[{\overline{Q}}^{-T}\overline{A}^{\left(2\right)}{}^T\overline{C}{\overline{A}}^{\left(2\right)}{\overline{Q}}^{-1}\right]\mathit{\hspace{1em}}\hfill & I_4\left[{\overline{Q}}^{-T}\overline{A}^{\left(2\right)}{}^T\overline{C}{\overline{A}}^{\left(2\right)}{\overline{Q}}^{-1}\right]\hfill \\ I_3\left[{\overline{Q}}^{-T}\overline{A}^{\left(2\right)}{}^T\overline{C}{\overline{A}}^{\left(1\right)}{\overline{Q}}^{-1}\right]\mathit{\hspace{1em}}\hfill & I_4\left[{\overline{Q}}^{-T}\overline{A}^{\left(2\right)}{}^T\overline{C}{\overline{A}}^{\left(2\right)}{\overline{Q}}^{-1}\right]\mathit{\hspace{1em}}\hfill & I_6\left[{\overline{Q}}^{-T}\overline{A}^{\left(2\right)}{}^T\overline{C}{\overline{A}}^{\left(2\right)}{\overline{Q}}^{-1}\right]\hfill \\ & & +I_7[{\overline{Q}}^{-T}\overline{A}^{\left(1\right)}{}^T\hfill \\ & & \overline{C}{\overline{A}}^{\left(1\right)}{\overline{Q}}^{-1}]\hfill \end{array}\right]$

(41)

3.4 Equation of Elementary Motion

By applying Lagrange’s equations Equation (7) to the system discretized by the DQFEM method showed above, we obtain the following system of differential equations:

$$[M^e]\{\ddot{q}\}+[K^e]\{\mathrm{q}\}=[0]$$

(42)

Where:

• $[M^e]$ and $[K^e]$ are respectively the elementary matrices of mass and stiffness.

• $\{\ddot{q}\}$ and $\{\mathrm{q}\}$ are respectively the global acceleration and displacement vectors suitable for DQFEM connectivity.

• Considering $[M]$ and $[K]$ respectively the total matrices after assembly of mass and stiffness, therefore the differential equation of motion becomes:

$$\left[\sum _e{M}^e\right]\{\ddot{q}\}+\left[\sum _e{K}^e\right]\{\mathrm{q}\}=[0]$$

(43)

The assembly of the global matrices is similar to that of the classic version of the finite elements method to ensure displacement and rotational compatibility at the nodes of adjacent elements. In the current application of DQFEM, boundary conditions are applied in the same way as the hp version of the finite element method.

4 Numerical Results and Discussion

In this part, numerical results are presented, in order to examine and discussing the free vibration behavior of FG-CNTs composite beams considering different types of boundary conditions, namely simply-supported (SS), clamped-clamped (CC), clamped-free (CF) and clamped-simply supported (CS). In the beginning, a validation studies are conducted to show exactness of the present model. Afterwards, a parametric study is made to inspect the influences of various physical and geometrical parameters. We have chosen in this investigation a Poly (methyl methacrylate), referred to as PMMA, is selected for the matrix with material properties $E_m=2.5$ GPa, ${\nu }_m=0.3$ and ${\rho }_m=1190$ kg/m${}^3$. In addition, the (10,10) armchair SWCNT is selected as the reinforcement and has the following properties (Shen 2009; Yas and Samadi 2012): $v_{cnt}=0.19$; ${\rho }_{cnt}=1400$ kg/m${}^3$; $E_{11}^{cnt}=600$ GPa; $E_{22}^{cnt}=10$ GPa and $G_{12}^{cnt}=17.2$ GPa. It is suitable to present the free vibration results in the dimensionless form as

$$\overline{\omega }=\omega A\sqrt{\frac{I_{00}}{A_{110}}}$$

where $A_{110}$ and $I_{00}$ are of beam constructed of pure matrix material, respectively.

4.1 Convergence Study of the Results

The inspection in Table 2 is for searching the appropriate number of terms (number of elements Ne and sampling points N) in the Differential Quadrature Finite Element Method (DQFEM) for the degree of convergence related to the 1st mode of natural frequencies ($\overline{\omega }$) of FG-CNTs composite beams. Various volume fraction of carbon nanotubes has been considered ($V_{cnt}^{*}$ $=$ 0.12, 0.17 and 0.28) immersed in the epoxy matrix with different shapes of reinforcement including (UD- -O, -X, -V). It is seen that, the convergent degree for the fundamental natural frequency is very fast using only Ne $=$ 4 and N $=$ 10. In addition, the choice of the combination of sampling number N with the number of elements Ne is free, which means that any combination can be chosen, more convergent the results can be obtained with an increase in the values of sampling points.

In Figures 3–3, a combination of the sampling number N with the number of elements Ne (Ne(i $+$ 1) $=$ Ne(i) $+$ 1 and Ne(i $+$ 1) $=$ Ne(i) $+$ 2) was used. At each iteration, a set of sampling points and two additional elements were added. We can see that the DQFEM method converge quickly, at $\mathit{\text{N}}=4$ and $\mathit{\text{Ne}}=1$ being the starting point of convergence, the results begin to converge from $\mathit{\text{N}}=6$ and $\mathit{\text{Ne}}=2$, which confirms the effectiveness of the proposed model, it gives a good result with a minimum of refinement.

To have good agreement with other works in the literature, and according to the convergence study of the DQFEM used in this work, we will choose a number of sampling $\mathit{\text{N}}=10$, and a number of element $\mathit{\text{Ne}}=2$ in the rest of the investigation.

Table 2 Convergence results of the first natural frequency of FG-CNTs beams with S-S boundary condition $(\mathit{\text{L}}/\mathit{\text{h}}=5)$

N	4			10			20
Ne	1			4			9
	$V_{cnt}^{*}$			$V_{cnt}^{*}$			$V_{cnt}^{*}$
Pattern	0.12	0.17	0.28	0.12	0.17	0.28	0.12	0.17	0.28
UD	1.6945	2.1650	2.4001	1.6430	2.0951	2.3314	1.6430	2.0951	2.3314
O	1.4867	1.9016	2.2030	1.4300	1.8219	2.1205	1.4300	1.8219	2.1205
X	1.7694	2.2516	2.4365	1.7231	2.1903	2.3759	1.7231	2.1903	2.3759
V	1.6948	2.1669	2.4145	1.6270	2.0730	2.3233	1.6270	2.0730	2.3233

Figure 2 Convergence of the first frequency of CNTs beams as a function of the sampling number of DQHFEM nodes N and mesh degree-Ne with (Ne(i) $=$ Ne(i) $+$ 1 and Ni $=$ Ni $+$ 2).

Figure 3 Convergence of the first frequency of CNTs beams as a function of the element number of DQHFEM nodes N and mesh degree-Ne with (Ne(i) $=$ Ne(i) $+$ 1 and Ni $=$ Ni $+$ 2).

4.2 Results Comparison

To evaluate precision of the obtained results predicted by the present method, the natural frequencies values of simply supported composite beam was compared with existing works in the literature in two stages. First, without considering CNTs reinforcement, a functionally graded (Ceramic-Metal) beam in which its material properties are $E_c=380$ GPa, ${\rho }_c=3960$ kg/m${}^3$, ${\nu }_c=0.3$ for ceramic and $E_m=70$ GPa, ${\rho }_m=2700$ kg/m${}^3$, ${\nu }_m=0.3$ for metal with the work provide Thai and Vo (2012) in Table 2, next stage is to consider a composite beam reinforced by carbon nanotubes (CNTs) in graded pattern with the results provided by (Wattanasakulpong and Ungbhakorn 2013; Tagrara et al. 2015).

Table 3 First three modes natural frequencies $\overline{\omega }$ of FG-beams. (N $=$ 10, Ne $=$ 2)

				p
L/h	Mode	Model		0	0.5	1	5	10
5	1		Thai and Vo (2012)	5.1527	4.4107	3.9904	3.4012	3.2816
			Present	5.1527	4.5328	4.2551	3.7055	3.4513
	2	TBT	Thai and Vo (2012)	17.8812	15.4588	14.0100	11.5431	11.0240
			Present	17.8812	15.3461	13.7723	11.2436	10.8418
	3		Thai and Vo (2012)	34.2097	29.8382	27.0979	21.7158	20.5561
			Present	30.2314	27.0488	24.8169	19.2826	17.9026
	1		Thai and Vo (2012)	5.1531	4.4110	3.9907	3.3998	3.2811
			Present	5.1531	4.5332	4.2554	3.7036	3.4509
	2	SBT	Thai and Vo (2012)	17.8868	15.4631	14.0138	11.5324	11.0216
			Present	17.8868	15.3501	13.7757	11.2345	10.8389
	3		Thai and Vo (2012)	34.2344	29.8569	27.1152	21.6943	20.5581
			Present	30.2314	27.0499	24.8182	19.2772	17.9040
20	1		Thai and Vo (2012)	5.4603	4.6511	4.2051	3.6485	3.5390
			Present	5.4603	4.7882	5.0022	4.0061	3.7385
	2	TBT	Thai and Vo (2012)	21.5732	18.3962	16.6344	14.3746	13.9263
			Present	21.5732	18.3864	16.6132	14.3429	13.9072
	3		Thai and Vo (2012)	47.5930	40.6526	36.7679	31.5780	30.5369
			Present	47.6040	40.7801	37.0412	31.8954	30.7086
	1		Thai and Vo (2012)	5.4603	4.6511	4.2051	3.6484	3.5389
			Present	5.4603	4.7882	4.5022	4.0060	3.7384
	2	SBT	Thai and Vo (2012)	21.5736	18.3965	16.6347	14.3728	13.9255
			Present	21.5736	18.3867	16.6134	14.3412	13.9065
	3		Thai and Vo (2012)	47.5950	40.6542	36.7692	31.5699	30.5337
			Present	47.6060	40.7818	37.0426	31.8871	30.7045

From the Table 3, it can be observed that our results match well with the results obtained by (Thai and Vo 2012) for various power law exponent (p), slenderness ratios (L/h) and mode number which confirm our model.

The accuracy of the presented free vibration results are compared with those of hyperbolic and sinusoidal order beam theories presented respectively by (Wattanasakulpong and Ungbhakorn 2013; Tagrara et al. 2015) for functionally graded reinforced CNT beam and the results are presented in Table 4. It is indicated that the actual model and solution procedure can accurately predict natural frequencies of FG-CNTs composite beams.

Table 4 Comparison of dimensionless fundamental frequencies for CNTs beams without elastic foundation ($\mathit{\text{L}}/\mathit{\text{h}}=15$, $V_{cnt}^{*}=0.12$)

Source		UD	O	X	V
HSDT	Source (*)	0.9745	0.7454	1.1151	0.8441
	Present	0.9745	0.7454	1.1152	0.9442
SSDT	Source (**)	0.9749	0.7446	1.1163	0.8443
	Source (*)	0.9745	0.7453	1.1152	0.8441
	Present	0.9745	0.7453	1.1152	0.9442
SSDT	Source (*)	0.9756	0.7440	1.1180	0.8448
	Present	0.9749	0.7446	1.1163	0.9444
*Wattanasakulpong, and Ungbhakorn (2013) **Tagrara et al. (2015).

4.3 Parametric Study

Figure 4 demonstrates the maximum values of non-dimensional frequencies against length to thickness ratios (L/h) of FG-CNTs reinforced beams under free vibration. The beams are supposed to be simply supported at both supports (S-S) and reinforced by $V_{cnt}^{*}=0.17$% of CNTs with various patterns distribution. It can be seen, that the highest natural frequency is detected in the case of the beam with FG-X reinforcement and then followed by free vibration of the beams with UD, FG-V, and FG-O of reinforcements, correspondingly. This is due to that the stiffness of the composite beam with –X distribution being much greater compared to those of others types.

To find out the significant influence of volume fractions of CNTs $V_{cnt}^{*}$ on dynamic response of FG composite beams with respect slenderness ratios (L/h) when FG is $-$X type and S-S board condition is depicted in Figure 5. Increasing volume fraction coefficient $V_{cnt}^{*}$ leads to increases in natural frequency of reinforced composite FG beam. Because, the Young’s modulus of becomes high and consequently the beam become more rigid.

Figure 4 Dimensionless fundamental frequencies $\overline{\omega }$ of CNT beam with various thickness ratios L/h ($V_{cnt}^{*}=0.17$, S-S).

Figure 5 Dimensionless fundamental frequencies of X-beam with various thickness ratios(S-S).

In Table 5, the obtained results of dimensionless natural frequency of the refined sinusoidal FG-CNTs beams considering various thickness ratios (L/h), different types distribution of CNTs beams and CNT volume fractions are presented for benchmark results. It is found that the FG-X beam is a stiffest one, and its frequency is high compared to other cases. With incrementing CNTs fraction for each pattern distribution into the beams, there is a significant increase of frequency. Moreover, it is seen that thickness ratio has more significant influence on frequency parameter when it takes a reduced value.

Table 5 Comparisons of fundamental frequencies $\overline{\omega }$ for CNTRC beams with different volume fraction $V_{cnt}^{*}$ (N $=$ 10, Ne $=$ 2, Reddy)

L/h	$V_{cnt}^{*}=0.12$				$V_{cnt}^{*}=0.17$				$V_{cnt}^{*}=0.28$
	UD	O	X	V	UD	O	X	V	UD	O	X	V
5	1.6430	1.4300	1.7231	1.6270	2.0951	1.8219	2.1903	2.0730	2.3314	2.1205	2.3759	2.3233
10	1.2582	1.0083	1.3891	1.2290	1.5689	1.2463	1.7364	1.5291	1.8259	1.5041	1.9637	1.7902
15	0.9745	0.7453	1.1152	0.9442	1.1983	0.9088	1.3761	1.1583	1.4361	1.1149	1.6114	1.3930
20	0.7807	0.5830	0.9123	0.7531	0.9529	0.7065	1.1168	0.9172	1.1605	0.8733	1.3368	1.1190

Table 6 Comparison of natural frequencies of different types of CNTRC beam for different boundary conditions $V_{cnt}^{*}=0.17$

Boundary	Type of CNT	L/h
Conditions	Distribution	5	10	15	20
S-S	UD-Beam	2.0951	1.5689	1.1983	0.9529
	O-Beam	1.8219	1.2463	0.9088	0.7065
	X-Beam	2.1903	1.7364	1.3761	1.1168
	V-Beam	2.0730	1.5291	1.1583	0.9172
C-C	UD-Beam	2.6514	2.2459	1.9559	1.7070
	O-Beam	2.3480	1.9546	1.6311	1.3715
	X-Beam	2.7356	2.3445	2.0941	1.8736
	V-Beam	2.6083	2.1423	1.8073	1.5370
C-S	UD-Beam	2.4035	1.9717	1.6247	1.3536
	O-Beam	2.1385	1.6552	1.2926	1.0404
	X-Beam	2.4800	2.1025	1.7956	1.5365
	V-Beam	2.3468	1.8499	1.4772	1.2061
C-F	UD-Beam	0.9718	0.6383	0.4594	0.3551
	O-Beam	0.7973	0.4842	0.3384	0.2584
	X-Beam	1.0504	0.7325	0.5416	0.4239
	V-Beam	0.8882	0.5508	0.3883	0.2975

The variation of natural frequencies of FG-CNTS reinforced composite beam based on (RSDBT) for different boundary conditions (S-S, C-S, C-C and C-F) are presented in Table 6, according to different distribution schemes of CNTs (UD, -X, -V, -O) and thickness ratios (L/h $=$ 5, 10, 15, 20) at ($V_{cnt}^{*}=0.17$).

By examining the results, it is found that the maximum non-dimensional frequency is obtained for FG-CNTs composite beam having C-C boundary condition due to stiffened edges, while the frequencies of clamped free (C-F) composite beam provides the lowest values. Another observation is that FG-X distribution has the highest natural frequencies compared to others patters distribution. In addition, we notice that as slenderness ratios (L/h) increase the natural frequencies reduce.

Figure 6 Influence of boundary conditions on the frequency parameter $\overline{\omega }$ of CNTs beam. (X-beam, $V_{cnt}^{*}=0.17$).

Figure 6 exhibits the effect of different boundary conditions on the dimensionless frequency of FG composite beam reinforced with CNTs with respect to thickness ratios, at $V_{cnt}^{*}=0.17$, and FG-X distribution. It is clear that natural frequency obtained for Clamped-Free boundary condition proceeds the lowest natural frequency, while it is intermediate for Simply-Simply then Clamped-Simply boundary conditions. However, the value of natural frequency is highest for the case of Clamped-Clamped end edges.

Figure 7 Variation of frequency parameter $\overline{\omega }$ S-S beam with various thickness ratios and different patterns CNTs distribution for $V_{cnt}^{*}=0.17$.

Figure 8 3D bar chart presentation of dimensionless fundamental frequencies (S-S) X-beam with various thickness ratios and volume fraction of carbon nanotubes.

Figure 7 displays a 3D bar chart form of nondimensional frequency ($\overline{\omega }$) versus various distribution shapes of reinforcement CNTs and slenderness ratios (L/h) for $V_{cnt}=0.17$. The results from this figure show that, the FG-X strengthened beams have the maximum natural frequency, which is the strongest composite beam followed by the UD-,V-and O-Beams, respectively. In addition, increasing the value of length to thickness ratio (L/h) leads to a decrement in the frequency parameter ($\overline{\omega }$), due to the reduction in area, and the stiffness of the beam decrease.

Figure 9 Vibration mode shapes of simply supported FG-X CNTs composite beam ($V_{cnt}^{*}=0.17$, L/h $=$ 5) using RSDBT.

The variation of the non-dimensional frequency for S-S reinforced CNTs beams is depicted in Figure 8 by bar chart through 3D representation, with respect to the different values of volume fraction of carbon nanotubes (V${}_{c\mathit{}n\mathit{}t}$ $=$ 0.12, 0,17 and 0,28) and length to thickness ratios (L/h) with FG-X dispersion pattern. It can be seen that the largest value of the CNTs leads to higher natural frequencies for each arrangement of reinforcement, due to of the enhancement of the local stiffness. It is also observed as before that the values of non-dimensional frequency reduce with increasing length to thickness parameter (L/h).

The first three mode shapes for the transverse displacements of the S-S FG-CNTs composite beam are illustrated in Figure 9 for FG-X pattern reinforcement, $V_{cnt}^{*}=0.17$, and $\mathit{\text{L}}/\mathit{\text{h}}=5$. One can notice that all vibration mode shapes exhibits dual coupled (shear-flexural) mode for present RSDBT. Thus, from this double coupled mode that highlights the effect of transverse shear stress on the vibration behavior of FG-CNTs beams. Therefore, this effect is important and must be taken into consideration in designing of thick composite beams.

5 Conclusions

Dynamic analysis of functionally graded carbon nanotube-reinforced composite beams is provided during this study based on refined higher order shear deformation beam theory, needless of the use shear correction factor. The functionally graded form, which describe the repartition of the carbon nanotubes across the thickness direction of the beam are considered, and their materials properties are estimated through extended rule of mixture. The general governing equations of motions have been extracted employing Lagrange’s principal and then solved via a robust numerical method, called differential quadrature finite element method for the first time in this area, with high rate of convergence speed and enhanced numerical stability. The exactitude of the present model in conjunction with the numerical resolution is verified by comparison with some existing results through a tabulated form. According to the obtained results from this investigation, the most interesting remarks can be summarized as follows:

• The present numerical tool based on (DQFEM) shows a fast convergence, rapid calculus in term of execution, and good numerical stability. The proposed model can be extended to study more complicated structures.

• The calculation times employing the present DQFEM based on refined shear deformation beam model is considerably faster compared to the classical finite element method and analytical solution.

• The convergence of the results can be controlled by increasing the number of samples and the number of elements.

• Overall, the presence of CNTs as being reinforcing elements with a graded form show an significant role in enhancing stiffness and strength of the composite beams.

• From the obtained results, it is shown that the FG-CNTs X pattern provides the higher non-dimensional natural frequency compared to other different types of CNTs beams, while the O-Beam has the smallest.

• It is revealed that, an increase in CNT volume fractions (V${}_{c\mathit{}n\mathit{}t}$) leads to an increment in natural frequencies of the reinforced composite beam.

• The change in length to thickness ratios (L/h) have significant on the maximum value of the natural frequency frequencies of FG-CNTs beams.

• The frequency results of the beams with C-C boundary condition are higher than those of C-S and S-S boundary conditions, respectively.

Acknowledgements

We acknowledge with grateful thanks the support by the laboratory of mechanical and material systems engineering, as well as the General Directorate of Scientific Research and Technological Development of the Ministry of Higher Education of Algeria.

References

Abdollahi, I., and Yas, M.H. 2020. Free vibration analysis of Timoshenko beams reinforced by BNNTs and a comparison with CNT-reinforced composite. SN Applied Sciences 2:645. doi:10.1007/s42452-020-2429-5.

Attia, M.A., and S.A. Mohamed. 2020. Nonlinear thermal buckling and postbuckling analysis of bidirectional functionally fraded tapered Microbeams Based on Reddy beam theory. Engineering with Computers. doi:10.1007/s00366-020-01080-1.

Azimi, M., Mirjavadi, S.S,. Navvab, S., A. M. S. Hamouda, and E. Davari. 2018. Vibration of rotating functionally graded Timoshenko nano-beams with nonlinear thermal distribution. Mechanics of Advanced Materials and Structures 25(6):467–480. doi:10.1080/15376494.2017.1285455.

Babaei, H., Kiani, Y., and Reza. Eslami, M. 2021. Vibrational behavior of thermally pre-/post-buckled FG-CNTRC beams on a nonlinear elastic foundation: a two-step perturbation technique. Acta Mechanica. doi:10.1007/s00707-021-03027-z.

Barati, M.R., and Shahverdi, H. 2020. Finite element forced vibration analysis of refined shear deformable nanocomposite graphene platelet-reinforced beams. Journal of the Brazilian Society of Mechanical Sciences and Engineering 42(33). doi:10.1007/s40430-019-2118-8.

Bekhadda, A., Bensaid, I., Cheikh, A., and Kerboua. B. 2019. Static Buckling and Vibration Analysis of Continuously Graded Ceramic-metal Beams Using a Refined Higher Order Shear Deformation Theory. Multidiscipline Modeling in Materials and Structures 15(6):1152–1169. doi:10.1108/MMMS-03-2019-0057/full/html.

Bensaid, I., and Saimi, A. 2022. Dynamic investigation of functionally graded porous beams resting on viscoelastic foundation using generalised differential quadrature method. Australian Journal of Mechanical Engineering. doi:10.1080/14484846.2021.2017115.

Ebrahimi, F., and Karimiasl, M. 2018. Nonlocal and surface effects on the buckling behavior of flexoelectric sandwich nanobeams. Mechanics of Advanced Materials and Structures 25(11):943–952. doi:10.1080/15376494.2017.1329468.

Eltaher, M.A., Mohamed, S.A., and Melaibari, A. 2020. Static Stability of unified composite Beams under Varying Axial Loads. Thin-Walled Structures 147:106488. doi:10.1016/j.tws.2019.106488.

Han, Y., and J. Elliott. 2007. Molecular dynamics simulations of the elastic properties of polymer / carbon nanotube composites. Computational Materials Science 39(2):315–323. doi:10.1016/j.commatsci.2006.06.011.

Iijima, S. 1991. Helical microtubules of graphitic carbon. Nature 354:56–58. doi:10.1038/354056a0.

Jam, J.E., and Kiani, Y. 2015. Low velocity impact response of functionally graded carbon nanotube reinforced composite beams in thermal environment. Composite Structures 132:35–43. doi:10.1016\%2Fj.compstruct.2015.04.045.

Karamanli, A., and Vo, T.P. 2021. Finite element model for carbon nanotube reinforced and graphene nanoplatelet-reinforced composite beams. Composite Structures 264:113739. doi:10.1016/j.compstruct.2021.113739.

Ke, L.L., Yang, J., and Kitipornchai, S. 2010. Nonlinear free vibration of functionally graded carbon nanotube-reinforced composite beams. Composite Structures 92(3):676–683. doi:10.1016/j.compstruct.2009.09.024.

Ke, L.L., Yang, J., and Kitipornchai, S. 2013. Dynamic stability of functionally graded carbon nanotube-reinforced composite beams. Mechanics of Advanced Materials and Structures 20(1):28–37. doi:10.1080/15376494.2011.581412.

Kiang. CH., M. Endo, P.M. Ajayan, G. Dresselhaus, and M.S. Dresselhaus. 1998. Size effects incarbon nanotubes. Physical Review Letters 81:1869–72. doi:10.1103/PhysRevLett.81.1869.

Kiani, Y. 2016. Shear buckling of FG-CNT reinforced composite plates using Chebyshev-Ritz method. Compos Part B Engineering 105:176–87. doi:10.1016/j.compositesb.2016.09.001.

Kiani, Y., and Mostafa, M. 2019. Nonlinear stability of sandwich beams with carbon nanotube reinforced faces on elastic foundation under thermal loading. Proc IMechE Part C: J Mechanical Engineering Science 233(5):1701–1712. doi:10.1177\%2F0954406218772613.

Lal, R., and Dangi, CH. 2019. Thermal vibrations of temperature-dependent functionally graded non-uniform Timoshenko nanobeam using nonlocal elasticity theory. Materials Research Express 6(7):075016. doi:10.1088/2053-1591/ab1332.

Lei, Z.X., Liew, K.M., and Yu, J.L. 2013. Free vibration analysis of functionally graded carbon nanotube-reinforced composite plates using the element-free kp-Ritz method in thermal environment. Composite Structures 106:128–138. doi:10.1016/j.compstruct.2013.06.003.

Lei, J., He, Y., Li, Z., Guo, S and Liu., D. 2019. Postbuckling analysis of bi-directional functionally graded imperfect beams based on a novel third-order shear deformation theory. Composite Structures. doi:10.1016/j.compstruct.2018.10.106.

Lin, F., and Xiang, Y. 2014. Vibration of carbon nanotube reinforced composite beams based on the first and third order beam theories. Applied Mathematical Modelling 38(15–16): 3741–3754. doi:10.1016/j.apm.2014.02.008.

Liu, C., Liu, B., Zhao, L., Xing, Y., Ma, C.H and Li, H. 2016. A differential quadrature hierarchical finite element method and its applications to vibration and bending of Mindlin plateswith curvilinear domains. International Journal for Numerical Methods in Engineering 109(2):174–197. doi:10.1002/nme.5277.

Saimi, A., A. Hadjoui, I. Bensaid, and A. Fellah. 2020. An Differential Quadrature Finite Element and the Differential Quadrature Hierarchical Finite Element Methods for the Dynamics Analysis of on Board Shaft. European Journal of Computational Mechanics 29 (4–6):303–344. doi:10.13052/ejcm1779-7179.29461.

Shen, H.S. 2009. Nonlinear bending of functionally graded carbon nanotube-reinforced composite plates in thermal environments. Composite Structures 91(1):9–19. doi:10.1016/j.compstruct.2009.04.026.

Shen, H.S., and Xiang, Y. 2013. Nonlinear analysis of nanotube reinforced composite beams resting on elastic foundations in thermal environments. Engineering Structures 56:698–708. doi:10.1016/j.engstruct.2013.06.002.

Tagrara, S., H., Benachour, A., Bouiadjra, M.B., and Tounsi, A. 2015. On bending, buckling and vibration responses of functionally graded carbon nanotube-reinforced composite beams. Steel and Composite Structures 19(5):1259–1277. doi:10.12989/scs.2015.19.5.1259.

Thai, H.T. and Vo, T.P. 2012. Bending and free vibration of functionally graded beams using various higher-order shear deformation beam theories. International Journal of Mechanical Sciences 62(1): 57–66. doi:10.1016/j.ijmecsci.2012.05.014.

Thostenson, E.T., Ren, Z.F., and Chou, T.W. 2001. Advances in the science and technology of carbon nanotubes and their composites: a review. Composites Science and Technology 61(13):1899–1912. doi:10.1016/S0266-3538(01)00094-X.

Vo-Duy, T., Ho-Huu, V., and Nguyen-Thoi, T. 2019. Free vibration analysis of laminated FG-CNT reinforced composite beams using finite element method. Frontiers of Structural and Civil Engineering 13:324–336. doi:10.1007/s11709-018-0466-6.

Vo, T., and Thai, H.T. Static behavior of composite beams using various refined shear deformation theories. Composite Structures 94(8):2513–2522. doi:10.1016/j.compstruct.2012.02.010.

Vo, T.P., Thai, H.T., and Aydogdu, M. 2017. Free vibration of axially loaded composite beams using a four-unknown shear and normal deformation theory. Composite Structures 178:406–414. doi:10.1016/j.compstruct.2017.07.022.

Vo, T.P., Thai, H.T., Nguyen, T.K., and F. Inam. 2014. Static and vibration analysis of functionally graded beams using refined shear deformation theory. Meccanica 49:155–168. doi:10.1007/s11012-013-9780-1.

Wang, Z.X., and Shen, H.S. 2011. Nonlinear vibration of nanotube-reinforced composite plates in thermal environments. Computational Materials Science 50(8):2319–2330. doi:10.1016/j.commatsci.2011.03.005.

Wattanasakulpong, N., and Ungbhakorn, V. 2013. Analytical solutions for bending, buckling and vibration responses of carbon nanotube-reinforced composite beams resting on elastic foundation. Computational Materials Science 71:201–208. doi:10.1016/j.commatsci.2013.01.028.

Wu, H., Kitipornchai, S., and Yang, J. 2015. Free vibration and buckling analysis of sandwich beams with functionally graded carbon nanotube-reinforced composite face sheets. International Journal of Structural Stability and Dynamics 15(7):1540011. doi:10.1142/S0219455415400118.

Xing, Y. F. and Liu, B. 2009. High-accuracy differential quadrature finite element method and its application to free vibrations of thin plate with curvilinear domain. International Journal for Numerical methods in engineering 80(13):1718–1742. doi:10.1002/nme.2685.

Yan, Y., Liu, B., Xing, Y., Carrera, E and Pagani, A. 2021. Free vibration analysis of variable stiffness composite laminated beams and plates by novel hierarchical differential quadrature finite elements. Composite Structures 274:114364. doi:10.1016/j.compstruct.2021.114364.

Yang, J., Ke, L.L., and Feng, C. 2015. Dynamic buckling of thermo-electro-mechanically loaded FG CNTRC beams. International Journal of Structural Stability and Dynamics 15(8):1540017. doi:10.1142/S0219455415400179.

Yarasca, J., Mantari, J.L., and Arciniega, R.A. 2016. Hermite–Lagrangian finite element formulation to study functionally graded sandwich beams. Composite Structures 140:567–581. doi:10.1016/j.compstruct.2016.01.015.

Yas, M.H., and Samadi, N. 2012. Free vibrations and buckling analysis of carbon nanotube-reinforced composite Timoshenko beams on elastic foundation. International Journal of Pressure Vessels and Piping 98:119–128. doi:10.1016\%2Fj.ijpvp.2012.07.012.

Biographies

Ihab Eddine Houalef, Ph.D student in Mechanical Engineering from Abou Beckr Belkaid University Tlemcen, Algeria. He is currently working in the level of the Mechanical engineering department at the same University. Mechanical and structural Engineering, Materials, Composite, Maintenance, Nanostructures and Dynamical Systems.

Ismail Bensaid received his B.Sc, M.Sc and Ph.D degrees in Mechanical Engineering from Abou Beckr Belkaid University Tlemcen, Algeria. He is currently working in the level of the Mechanical engineering department at the same University. Dr. Bensaid does research in Mechanical and structural Engineering, Materials, Composite, Maintenance, Nanostructures and Dynamical Systems. He, as an author/co-author, has published more than 18 articles in various journals.

Ahmed Saimi obtained his Ph.D in Mechanics of Materials and Structures from the University of Tlemcen, Algeria, in 2017. He is currently a Senior Lecturer at the National High School of Hydraulics Blida, Algeria. A researcher member of Mechanical Systems and Structural Engineering Laboratory, IS2M/UABT. His research interests are: Finite element methods, Structural vibration, Structural dynamics, Dynamics of rotors, Dynamics of rotating machines, computational mechanics, FG materials, Composite materials.

Abdelmadjid Cheikh obtained his Ph.D in mechanical engineering from the University of Tlemcen, Algeria. He is currently a professor at the University of Tlemcen, Algeria (UABT). Research Director in Mechanical Systems and Structural Engineering Laboratory, IS2M/UABT. His research interests are: Materials Engineering, Structural Engineering, Mechanical Engineering, Structural Analysis, Finite elements Modeling, Structural Dynamics, Simulation, Dynamic Analysis, Modal Analysis, Structural Vibration, Vibration Analysis, mechanical fabrication, tolerancing, CAO, DAO.

European Journal of Computational Mechanics, Vol. 31_4, 505–538.
doi: 10.13052/ejcm2642-2085.3142
© 2023 River Publishers