Vibration Analysis of Fluid Conveying Carbon Nanotubes Based on Nonlocal Timoshenko Beam Theory by Spectral Element Method

In this work, we applied the spectral element method (SEM) to analyze the dynamic characteristics of fluid conveying single-walled carbon nanotubes (SWCNTs). First, the dynamic equations for fluid conveying SWCNTs were deduced based on the nonlocal Timoshenko beam theory. Then, the spectral element formulation was established for a free/forced vibration analysis of fluid conveying SWCNTs by introducing discrete Fourier transform. Furthermore, the proposed method was validated using several comparison examples. Finally, the natural frequencies and dynamic responses of a simply-supported fluid conveying SWCNTs were calculated by the SEM, considering different internal fluid velocities and small-scale parameters (SSPs). The effects of fluid velocity and SSPs on the dynamic characteristics of SWCNTs conveying fluid were revealed by the numerical results. Compared with other methods, the SEM shows high accuracy and efficiency.


Introduction
Carbon nanotubes (CNTs) have been extensively used in numerous areas due to their excellent mechanical properties [1,2]. A deep understanding of the mechanical behavior of CNTs leads to wider applications [3,4]. In general, a molecular dynamics simulation and experimental methods are usually adopted to predict the mechanical behavior of a CNT [5][6][7]. Since the experiments at the nanoscale level are difficult to conduct [8] and have high computational cost, the molecular dynamics simulation method is more often used to deal with small-size atomic systems. Hence, the continuum-based solid mechanics models or methods [8][9][10][11] bring vitality for predicting the mechanical properties of a large CNT. As an important mechanical property, the vibration characteristics of a CNT have become another hot topic for research in recent years. For example, Cai et al. [12,13] discussed the dynamic responses of a CNT-based nanobeam, which was considered a model of nanobalance. For slight vibrations, the deformation is linearly elastic and reversible. Therefore, elastic beam models can be adopted to analyze the dynamic properties of a CNT [14]. The conventional Euler-Bernoulli and Timoshenko beam theories, however, cannot be used directly to predict the vibration characters of CNTs because of the scale effect; this defect was overcome by the nonlocal elasticity theory proposed by Eringen et al. [15,16]. Based on nonlocal beam theory, the static bending, buckling, vibration, and wave properties of single-/double-walled CNTs have been studied systematically [17][18][19][20][21]. In recent years, the nonlocal mechanics of elastic nanobeams have undergone rapid progress, with a lot of newly-developed theories, e.g., strain-driven and stress-driven nonlocal integral elasticity theory [22][23][24], two-phase integral elasticity theory [25], nonlocal strain gradient elasticity theory [26][27][28], and modified nonlocal strain gradient elasticity theory [29,30].

Dynamic Equations of SWCNT Conveying Fluid
The fluid conveying SWCNT with a free body diagram is shown in Figure 1. By assuming incompressible and non-viscous uniform flow, and considering the coupling effects (the centrifugal and Coriolis forces of fluid) between internal fluid and tube, the transverse vibration equation can be written as follows [46]: where V is the velocity of the internal flow and y is the transverse deflection of the tube, i.e., the function of axial coordinate x and time t.
where V is the velocity of the internal flow and y is the transverse deflection of the tube, i.e., the function of axial coordinate x and time t.
According to the nonlocal elastic theory, the dynamic equation of micro-Timoshenko beam takes the following form [19]: where E is Young's modulus, G is the shear modulus, I is the moment of inertia of the beam cross section about its neutral axis, k0 is the shear coefficient in Timoshenko beam theory, μ is the Poisson ratio, Ap is the area of the beam cross section, Jp is the mass moment of inertia, mp is the mass of tube per unit length, mf is the mass of fluid per unit length, and e0a is the scale coefficient to reflect a smallscale effect (Note: e0 is the material constant obtained by experiments or molecular dynamics, and a is the internal characteristic length, e.g., for CNT, the C-C bond length is usually chosen as the internal characteristic length [8,19]), Q is the shear force, M is the bending moment, and φ is the rotation angle of tube's cross section caused by bending deformation.
According to Equations (1) and (2), the following partial differential equation can be obtained after eliminating Q and M: By introducing the following non-dimensional parameters, the aforementioned equations will be relevantly simplified: According to the nonlocal elastic theory, the dynamic equation of micro-Timoshenko beam takes the following form [19]: where E is Young's modulus, G is the shear modulus, I is the moment of inertia of the beam cross section about its neutral axis, k 0 is the shear coefficient in Timoshenko beam theory, µ is the Poisson ratio, A p is the area of the beam cross section, J p is the mass moment of inertia, m p is the mass of tube per unit length, m f is the mass of fluid per unit length, and e 0 a is the scale coefficient to reflect a small-scale effect (Note: e 0 is the material constant obtained by experiments or molecular dynamics, and a is the internal characteristic length, e.g., for CNT, the C-C bond length is usually chosen as the internal characteristic length [8,19]), Q is the shear force, M is the bending moment, and ϕ is the rotation angle of tube's cross section caused by bending deformation.
According to Equations (1) and (2), the following partial differential equation can be obtained after eliminating Q and M: By introducing the following non-dimensional parameters, the aforementioned equations will be relevantly simplified: where L is the length of the CNT. Substituting the aforementioned terms in Equation (5) into Equations (3) and (4), one can obtain the dimensionless equations as follows: Nanomaterials 2019, 9, 1780 4 of 16 Taking the pairs of the discrete Fourier transform (Equation (8)) of η and ϕ back into Equations (6) and (7), the partial differential equations in the time domain can be transformed into ordinary differential equations in the frequency domain: whereη n is the discrete Fourier transform pair of η r , ω is the dimensionless form of frequency Ω, and the subscript in ω n presents the discrete dimensionless frequency. N denotes the number of samples adopted in the discrete Fourier transform, which is a kind of series expansion method [47]. In general, a more efficient fast Fourier transform (FFT) algorithm is used during the application. Accordingly, N needs to satisfy N = 2 n and n is a positive integer, τ r is the discrete dimensionless time, and i is an imaginary unit. Substituting Equation (8) into Equations (6) and (7), one can get the ordinary differential equations with respect toη n andφ n . Since these equations must be satisfied at every discrete frequency point, the subscript n can be omitted, and Equations (6) and (7) become The solution of the ordinary differential Equations (9) and (10) with constant coefficients holds the following exponential form:η (ξ) = Ce ikξ ,φ(ξ) = De ikξ (11) where k is the dimensionless wave number, and C, D are undetermined coefficients. Substituting Equation (11) into Equations (9) and (10), the following matrix equation can be obtained: The dispersion relation of elastic wave propagation in SWCNT conveying fluid can be obtained using the nontrivial solution condition of Equation (12), i.e., The four roots of polynomial Equation (13) represent four dimensionless wave numbers. In general, half of them are real roots representing the propagating wave motion, and the other half are imaginary roots corresponding to evanescent wave motion. For one-dimensional wave propagation, one may use k 1 , k 2 denoting a left-traveling wave, and k 3 , k 4 a right-traveling wave. Due to the four wave numbers, Equation (10) can also be reformulated as follows: where D j = λ j C j , and the coefficient λ j can be determined by Equation (12), i.e., Similarly, both the dimensionless shear force and bending moment can be expressed aŝ with Now, all the related parameters in the frequency-domain formulations are obtained. In the following section, the spectral formulation of a finite length SWCNT conveying fluid will be established.

The Spectral Formulations
The one-dimensional wave motion in SWCNT ( Figure 2) includes left-and right-traveling waves w l and w r . For a SWCNT element with length l, Equation (14) holds the following form if node 1 was chosen as origin [45]: Nanomaterials 2019, 9, x FOR PEER REVIEW 6 of 18

The Spectral Formulations
The one-dimensional wave motion in SWCNT ( Figure 2) includes left-and right-traveling waves l w and r w . For a SWCNT element with length l, Equation (14) holds the following form if node 1 was chosen as origin [45]: The nodal displacement (Figure 3a) of the element can be practically represented as.
The nodal displacement (Figure 3a) of the element can be practically represented as.

The Spectral Formulations
The one-dimensional wave motion in SWCNT ( Figure 2) includes left-and right-traveling waves l w and r w . For a SWCNT element with length l, Equation (14) holds the following form if node 1 was chosen as origin [45]: The nodal displacement (Figure 3a) of the element can be practically represented as. Putting a specific position parameter into Equations (18) and (19) can be formulated in a matrix form, i.e., Putting a specific position parameter into Equations (18) and (19) can be formulated in a matrix form, i.e., Substituting the coefficient column in Equation (20) into Equation (18) yields another expression of nodal displacements as follows: In fact, Equation (22) presents the shape function of the spectral element. Similarly, the nodal force shown in Figure 3b reads In the same way, the nodal force can be written in a matrix form as follows: Based on Equations (20) and (24), the relationship between the nodal force and nodal displacement can be obtained, i.e., where matrix S (e) = BA −1 is called the spectral element matrix. The natural frequency of the system can be obtained when Equation (26) has nontrivial solution, i.e., det S (e) = 0.
For a simply supported SWCNT, both ends are chosen as the element nodes, and the boundary conditions are as follows: Substituting the aforementioned boundary conditions into Equation (26) and picking the related four elements in the matrix, one can get the following homogeneous equation: In Equation (28)

Comparison Example
To validate the SEM, we calculated the first five natural frequency parameters √ ω of SWCNT with two different boundary conditions and compared the results with those in reference [19]. The parameters of SWCNT are listed in Table 1, and the results of simply supported SWCNTs with different SSPs α are listed in Table 2.   There are several methods for computing the shear coefficient [48,49]. In this paper, the value of k 0 can be obtained by k 0 = 6(1+µ)(1+m) 2 (7+6µ)(1+m) 2 +(20+12µ)m 2 for a Timoshenko beam with a hollow circle cross section [50,51].
It can be seen in Table 2 that the results in this paper are the same as those in reference [19]. Therefore, the SEM is available to predict the natural frequencies of a SWCNT in free vibration.  We also computed the first five natural frequency parameters  of a cantilevered SWCNT with four different SSPs α; the comparison results are shown in Table 3. The results by the SEM and in reference [19] are identical.   We also computed the first five natural frequency parameters √ ω of a cantilevered SWCNT with four different SSPs α; the comparison results are shown in Table 3. The results by the SEM and in reference [19] are identical. for the cantilevered SWCNT with different dimensionless SSPs α are shown in Figure 5a,b. For the cantilevered SWCNT, the effect of α on the natural frequencies is slightly different from that on a simply supported SWCNT. It can be seen in Figure 5 that all the natural frequencies decrease with increasing α, except the first-order natural frequency, which actually increases. four different SSPs α; the comparison results are shown in Table 3. The results by the SEM and in reference [19] are identical.  Figure 5 (a,b). For the cantilevered SWCNT, the effect of α on the natural frequencies is slightly different from that on a simply supported SWCNT. It can be seen in Figure 5 that all the natural frequencies decrease with increasing α, except the first-order natural frequency, which actually increases.

Free Vibration of a SWCNT Conveying Fluid
In this example, we analyzed the free vibrations of a simply supported SWCNT conveying fluid with the same material parameters as those in Table 1. The density of the internal fluid was 1.0 g/cm 3 . For the purpose of investigating the effects of α and fluid velocity u on the natural frequencies, we calculated the first five natural frequencies of SWCNT conveying fluid in two cases, i.e., constant fluid velocity u = 0.4 with different SSPs α and constant α = 0.2 with different fluid velocities. The results of the former case are listed in Table 4, and those of the latter case are listed in Table 5. As shown in Figure 6, the effect of α on the SWCNT conveying fluid is the same as that on the hollow SWCNT. All the natural frequencies decrease with increasing α.
For the purpose of investigating the effects of α and fluid velocity u on the natural frequencies, we calculated the first five natural frequencies of SWCNT conveying fluid in two cases, i.e., constant fluid velocity u = 0.4 with different SSPs α and constant α = 0.2 with different fluid velocities. The results of the former case are listed in Table 4, and those of the latter case are listed in Table 5. As shown in Figure 6, the effect of α on the SWCNT conveying fluid is the same as that on the hollow SWCNT. All the natural frequencies decrease with increasing α.   The effect of the internal fluid velocity on the natural frequencies of SWCNT conveying fluid, as shown in Figure 7, is the same as that of a macropipe conveying fluid. The flowing fluid weakens the tube effective stiffness and leads to a decrease in the natural frequencies. When the first natural frequency drops to zero, divergence (buckling) instability occurs; the corresponding velocity is named critical velocity. The effect of the internal fluid velocity on the natural frequencies of SWCNT conveying fluid, as shown in Figure 7, is the same as that of a macropipe conveying fluid. The flowing fluid weakens the tube effective stiffness and leads to a decrease in the natural frequencies. When the first natural frequency drops to zero, divergence (buckling) instability occurs; the corresponding velocity is named critical velocity. Moreover, the critical velocity of the simply supported fluid conveying SWCNTs was computed with the mass ratio β = 0.64, α = 0.3 to compare with that in reference [37]. The variation of the fundamental frequency with fluid velocity is shown in Figure 8. Herein, the critical velocity uc corresponding to zero fundamental frequency is uc = 2.254, which is slightly lower than that predicted by Wang [37] based on the nonlocal Euler-Bernoulli beam model with uc = 2.286. Moreover, the critical velocity of the simply supported fluid conveying SWCNTs was computed with the mass ratio β = 0.64, α = 0.3 to compare with that in reference [37]. The variation of the fundamental frequency with fluid velocity is shown in Figure 8. Herein, the critical velocity u c corresponding to zero fundamental frequency is u c = 2.254, which is slightly lower than that predicted by Wang [37] based on the nonlocal Euler-Bernoulli beam model with u c = 2.286. Moreover, the critical velocity of the simply supported fluid conveying SWCNTs was computed with the mass ratio β = 0.64, α = 0.3 to compare with that in reference [37]. The variation of the fundamental frequency with fluid velocity is shown in Figure 8. Herein, the critical velocity uc corresponding to zero fundamental frequency is uc = 2.254, which is slightly lower than that predicted by Wang [37] based on the nonlocal Euler-Bernoulli beam model with uc = 2.286.

Spectral Formulations
A simply supported SWCNT conveying fluid, as shown in Figure 9, was subjected to a point load in the middle part. w denotes the wave motion, while the subscripts l and r denote the left and right part of the load position or wave motion, respectively.

Spectral Formulations
A simply supported SWCNT conveying fluid, as shown in Figure 9, was subjected to a point load in the middle part. w denotes the wave motion, while the subscripts l and r denote the left and right part of the load position or wave motion, respectively. The SWCNT is divided into two elements when an external point load exists at the middle node p. Accordingly, the spectral matrix in Equation (26) is rewritten in the following global form: The nodal force and displacement column vectors are changed into the following form: For simplicity, the sequence of nodal force and displacement column vectors were rearranged as the following formulas [52]: The SWCNT is divided into two elements when an external point load exists at the middle node p. Accordingly, the spectral matrix in Equation (26) is rewritten in the following global form: The nodal force and displacement column vectors are changed into the following form: For simplicity, the sequence of nodal force and displacement column vectors were rearranged as the following formulas [52]: It is obvious that the relationship between the original and the current nodal force and displacement column vectors can be expressed in the following form by introducing the reordered matrix R [52]: The simply supported boundary conditions can be modified in the following matric forms by introducing the boundary matrices T b and T f [52]: Now, the reordered nodal force and displacement vectors in Equation (22) can be rewritten as Substituting Equations (33) and (35) into Equation (26) and multiplying both sides by R T , the following results can be obtained: With S R = R T S (e) R .

Multiplying both sides of Equation (36) byT
T b and considering the result T T b T f = 0 2×2 [52], one can obtain the following result: The unknown nodal displacement vector can be resolved from Equation (37) as follows [52]: Once the nodal displacement is obtained, the displacement at any location can be obtained according to the shape function expression of Equation (22): with w 1 = 0φ 1ηpφp T , w 2 = η pφp 0φ 2 T . Equation (39) is the dynamic response of the SWCNT conveying fluid in the frequency domain, and that in the time domain can be quickly obtained through the inverse discrete Fourier transform of Equation (39).

Comparison Example
We also provided an example to demonstrate the validity of the SEM in dynamic response prediction. Up to now, no example of microbeams has been given; as such, the macro-Timoshenko beam was chosen here for comparison [53]. The beam with a rectangular cross-section is simply supported at both ends, and the point harmonic load is exerted at the center of the beam. The material parameters are listed in Table 6.
The steady harmonic load is F(x, t) = 10 6 δ(x − L/2) sin(1000t), and the transverse displacement time history curves at x = L/2 are shown in Figure 10. Obviously, the results by the SEM are identical to those by the boundary element method (BEM) [53].

Comparison Example
We also provided an example to demonstrate the validity of the SEM in dynamic response prediction. Up to now, no example of microbeams has been given; as such, the macro-Timoshenko beam was chosen here for comparison [53]. The beam with a rectangular cross-section is simply supported at both ends, and the point harmonic load is exerted at the center of the beam. The material parameters are listed in Table 6.
The steady In reference [53], the time step was chosen as 6 10 5    t . However, the results with the same accuracy can be obtained using 4 10 1    t in this work. Therefore, the SEM not only holds a high level of accuracy, but also takes low computational cost in dynamic analysis.

Example of SWCNT Conveying Fluid
We finally studied the dynamic response of a simply supported SWCNT conveying fluid under non-dimensional harmonic load  In reference [53], the time step was chosen as ∆t = 5 × 10 −6 . However, the results with the same accuracy can be obtained using ∆t = 1 × 10 −4 in this work. Therefore, the SEM not only holds a high level of accuracy, but also takes low computational cost in dynamic analysis.

Example of SWCNT Conveying Fluid
We finally studied the dynamic response of a simply supported SWCNT conveying fluid under non-dimensional harmonic load F(ξ, τ) = δ(ξ − 1/2) sin(πτ) exerted on the midpoint of the tube.
The same parameters of the same tube and internal fluid as those in Section 3.3 were used. Two cases with u = 0.4 and α = 0.2 are considered.
When u = 0.4, the dimensionless transverse displacement at the center of the tube is drawn in Figure 11 with different α values. As seen in Figure 11, the dynamic responses become larger as the SSP α increases. This effect of α on the SWCNT vibrations is the same as that in Section 3.3 for the free vibration situation. Such a phenomenon reflects the fact that the SWCNT will become more and more flexible if α increases.
When α = 0.2, the time history of the dimensionless displacements is shown in Figure 12 for SWCNT conveying fluid with different fluid velocities. In Figure 12, one can find that the internal fluid velocity has a significant influence on the dynamic response of SWCNT. This phenomenon can also be observed in the vibration of a macropipe conveying fluid. This influence, however, will not be as significant as that in a macropipe. The reason for this is that the mass of the fluid in SWCNT is lower, and the CNT has a relatively higher bending stiffness. So the instability of the SWCNT conveying fluid occurs only when the internal flow has an extremely high velocity.

Concluding Remarks
The SEM for the analysis of flow-induced vibration of a SWCNT conveying fluid was developed based on the nonlocal Timoshenko beam theory. The spectral formulations for free/forced vibration of SWCNT conveying fluid were established. Numerical examples were given to show the validity of the modified method. Some conclusions were drawn, as follows.  As seen in Figure 11, the dynamic responses become larger as the SSP α increases. This effect of α on the SWCNT vibrations is the same as that in Section 3.3 for the free vibration situation. Such a phenomenon reflects the fact that the SWCNT will become more and more flexible if α increases.
When α = 0.2, the time history of the dimensionless displacements is shown in Figure 12 for SWCNT conveying fluid with different fluid velocities. As seen in Figure 11, the dynamic responses become larger as the SSP α increases. This effect of α on the SWCNT vibrations is the same as that in Section 3.3 for the free vibration situation. Such a phenomenon reflects the fact that the SWCNT w In Figure 12, one can find that the internal fluid velocity has a significant influence on the dynamic response of SWCNT. This phenomenon can also be observed in the vibration of a macropipe conveying fluid. This influence, however, will not be as significant as that in a macropipe. The reason for this is that the mass of the fluid in SWCNT is lower, and the CNT has a relatively higher bending stiffness. So the instability of the SWCNT conveying fluid occurs only when the internal flow has an extremely high velocity.

Concluding Remarks
The SEM for the analysis of flow-induced vibration of a SWCNT conveying fluid was developed based on the nonlocal Timoshenko beam theory. The spectral formulations for free/forced vibration of SWCNT conveying fluid were established. Numerical examples were given to show the validity of the modified method. Some conclusions were drawn, as follows.  In Figure 12, one can find that the internal fluid velocity has a significant influence on the dynamic response of SWCNT. This phenomenon can also be observed in the vibration of a macropipe conveying fluid. This influence, however, will not be as significant as that in a macropipe. The reason for this is that the mass of the fluid in SWCNT is lower, and the CNT has a relatively higher bending stiffness. So the instability of the SWCNT conveying fluid occurs only when the internal flow has an extremely high velocity.

Concluding Remarks
The SEM for the analysis of flow-induced vibration of a SWCNT conveying fluid was developed based on the nonlocal Timoshenko beam theory. The spectral formulations for free/forced vibration of SWCNT conveying fluid were established. Numerical examples were given to show the validity of the modified method. Some conclusions were drawn, as follows.
First, by calculating the first five natural frequencies for a SWCNT and the dynamic response for a macro-Timoshenko beam, the SEM is proved to be effective to predicts the dynamic characteristics of a SWCNT.
Second, the effect of α on the SWCNT conveying fluid is the same as that on the hollow SWCNT. All the natural frequencies decrease with increasing α.
Third, the effect of the internal fluid velocity on the natural frequencies of SWCNT conveying fluid is the same as that of a macropipe conveying fluid. The flowing fluid weakens the tube effective stiffness and leads to a decrease in the natural frequencies. When the flow velocity approaches the critical value, the first natural frequency drops to zero, and divergence (buckling) instability occurs.
Finally, compared with the BEM, the SEM can provide results with the same level of accuracy with a much larger time step, e.g., the time step in the present method is~500 times higher than that in the BEM.