*Article* **Free Vibration of Thin-Walled Composite Shell Structures Reinforced with Uniform and Linear Carbon Nanotubes: Effect of the Elastic Foundation and Nonlinearity**

**Avey Mahmure <sup>1</sup> , Francesco Tornabene 2,\* , Rossana Dimitri <sup>2</sup> and Nuri Kuruoglu <sup>3</sup>**


**Abstract:** In this work, we discuss the free vibration behavior of thin-walled composite shell structures reinforced with carbon nanotubes (CNTs) in a nonlinear setting and resting on a Winkler– Pasternak Foundation (WPF). The theoretical model and the differential equations associated with the problem account for different distributions of CNTs (with uniform or nonuniform linear patterns), together with the presence of an elastic foundation, and von-Karman type nonlinearities. The basic equations of the problem are solved by using the Galerkin and Grigolyuk methods, in order to determine the frequencies associated with linear and nonlinear free vibrations. The reliability of the proposed methodology is verified against further predictions from the literature. Then, we examine the model for the sensitivity of the vibration response to different input parameters, such as the mechanical properties of the soil, or the nonlinearities and distributions of the reinforcing CNT phase, as useful for design purposes and benchmark solutions for more complicated computational studies on the topic.

**Keywords:** CNT; elastic foundations; nonlinear free vibration; nonlinear frequency; shallow shell structures

#### **1. Introduction**

The fast development of nanotechnology in recent years has encouraged the production of nanotubes, increasing their application in many engineering areas. The CNTs produced for the first time by Iijima in 1993, are increasingly used in various industries and commonly proposed as novel material due to their great potential [1,2]. One of the most important application areas of CNTs stems from their large use as reinforcement phase in traditional composites and polymers. The mechanical, thermal and electrical properties of composites reinforced with CNTs are significantly improved compared to more classical composites, along with an increased level of strength in their structural application [3–5]. For such reasons, CNTs are used in some areas of the defense industry, especially in rocket, aerospace and aviation industries, where high-precision computations are required [6–8]. Among various problems is the linear and nonlinear vibration behavior of composite shell structures involving the presence of different distributions of CNTs. Composite shells, indeed, can include uniform or nonuniform patterns of CNTs, depending on the desired mechanical properties of the structures [9–24]. In this framework, a pioneering work on the nonlinear vibrations of composite shell structures was represented by [9], which considered a linear distribution of CNTs within the material. Following this work, some linear and nonlinear free vibration problems were proposed in [10–17] and [18–24], respectively, for

**Citation:** Mahmure, A.; Tornabene, F.; Dimitri, R.; Kuruoglu, N. Free Vibration of Thin-Walled Composite Shell Structures Reinforced with Uniform and Linear Carbon Nanotubes: Effect of the Elastic Foundation and Nonlinearity. *Nanomaterials* **2021**, *11*, 2090. https://doi.org/10.3390/ nano11082090

Academic Editor: Yang Tse Cheng

Received: 21 July 2021 Accepted: 15 August 2021 Published: 17 August 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

unconstrained shallow shells and panels reinforced by CNTs, while proposing different numerical methods to solve the related problems.

The technological evolution of artificial materials and their manufacturing has expanded the application areas for such materials, improving the interest towards even more complicated and coupled problems, as well as the possible interactions of a structural member with its surrounding medium. Composite CNT-based shell structures resting on elastic foundations can be found in different civil and mechanical engineering applications, in nuclear power plants, etc. Among different possibilities to model an elastic foundation, the Pasternak and Winkler models are two common ways of describing its mechanical behavior, based on a different number of input parameters [25,26]. When modeling the vibrations of structural members containing CNTs, it is important to study the effect of the reinforcement phase and elastic foundations on the frequency–amplitude relationships. Up to date, most works from the literature have been devoted to the solution of linear vibration problems, by means of different numerical techniques [27–35]. More specifically, Tornabene et al. [27] examined the influence of Winkler–Pasternak foundations on the static and dynamic analysis of laminated double-curved shells and panels using the differential quadrature method. The same numerical approach was successfully proposed in [28] to study the vibration response of functionally graded carbon nanotube reinforced composite (FG-CNTRC) spherical shells on an elastic foundation. Zhang and Liew [29] applied an element-free approach to study the large deflection response of FG-CNTRC plates. Dinh and Nguyen [30] applied a fourth-order Runge–Kutta method and Galerkin method to solve the dynamic and vibration problem of FG-CNTRC truncated conical shells on elastic foundations. Shen and He [31] performed a large amplitude vibration analysis of FG-CNTRC double-curved panels on elastic foundation by applying a two-step perturbation approach, as also implemented in [32] to analyze the large amplitude vibration of FG shallow arches on a nonlinear elastic foundation. A further linear formulation was proposed by Sobhy and Zenkour [33] to study the vibrations of FG graphene platelet reinforced composite double-curved shallow shells on an elastic foundation; Sofiyev et al. [34,35] investigated the stability of CNTRC conical shells resting on an elastic foundation under hydrostatic pressure and combined loads in different settings.

Despite the considerable attention paid by the scientific literature to the linear vibration of shell structures, the nonlinear vibrations of CNT shallow shells resting on elastic foundations have not been adequately investigated. In this context, this paper aims to study the nonlinear free vibration behavior of thin-walled shell structures reinforced with CNTs and resting on an elastic Winkler- or Pasternak-type foundation, while proposing a Grigolyuk method to handle the problem. The organization of the rest of the paper is as follows: Section 2 recalls the basic theoretical aspects for both the shell-foundation interaction and nonlinear structural problem. Section 3 illustrates the analytical methodology applied to solve the problem, whose numerical investigation is presented and discussed in Section 4, while Section 5 closes the work with main comments and remarks.

#### **2. Theoretical Formulation**

#### *2.1. Description of Shell-Foundation Interaction Model*

Let us consider a composite spherical and hyperbolic paraboloidal (hypar) shallow shell reinforced with CNTs with length *a*, width *b*, thickness *h* and curvature radii *R*<sup>1</sup> and *R*2, respectively (see Figure 1a,b). The Cartesian coordinate system (*x*1, *x*2, *x*3) is here assumed to define the shell geometry in its length, width and thickness direction, respectively. As also shown in Figure 1, both the spherical and hypar shallow shells are immersed in an elastic WPF, here modeled as follows [25,26]:

$$K(w) = k\_w w - k\_p \left(\frac{\partial^2 w}{\partial x\_1^2} + \frac{\partial^2 w}{\partial x\_2^2}\right) \tag{1}$$

where *kw* (in Pa/m) is the Winkler spring stiffness and *kP* (in *Pa* · *m*) refers to the shear layer stiffness. When *kP* = 0, the foundation reverts to a Winkler-type elastic foundation (WF). The FG-CNTRC shell structures feature the following properties [9]

$$\begin{aligned} \mathbf{Y}\_{11}^{\overline{\mathbf{T}}\_{3}} &= \eta\_{1} V\_{\text{CN}}^{\overline{\mathbf{T}}\_{3}} Y\_{11}^{\text{CN}} + V\_{\text{W}} E^{\text{m}}, \quad \frac{\eta\_{2}}{Y\_{22}^{\overline{\mathbf{T}}\_{3}}} = \frac{V\_{\text{CN}}^{\overline{\mathbf{T}}\_{3}}}{Y\_{22}^{\text{CN}}} + \frac{V\_{\text{m}}}{E^{\text{m}}}, \quad \frac{\eta\_{3}}{Y\_{12}^{\overline{\mathbf{T}}\_{3}}} = \frac{V\_{\text{CN}}^{\overline{\mathbf{T}}\_{3}}}{Y\_{12}^{\text{CN}}} + \frac{V\_{\text{m}}}{Y^{\text{m}}}, \\\\ \boldsymbol{\nu}\_{12} &= V\_{\text{CN}}^{\ast} \boldsymbol{\nu}\_{12}^{\text{CN}} + V\_{\text{m}} \boldsymbol{\nu}^{\text{m}}, \; \rho\_{1}^{\overline{\mathbf{T}}\_{3}} = V\_{\text{CN}}^{\ast} \rho^{\text{CN}} + V\_{\text{m}} \rho^{\text{m}}, \; \overline{\mathbf{x}}\_{3} = \boldsymbol{x}\_{3}/\hbar \end{aligned} \tag{2}$$

where the elastic properties for CNTs and matrix denoted as *YCN ij* (*i*, *<sup>j</sup>* = 1, 2), and *<sup>Y</sup>m*,*Gm*, respectively; *<sup>η</sup>j*(*<sup>j</sup>* = 1, 2, 3) refers to the efficiency parameters for CNTs; *<sup>V</sup>x*<sup>3</sup> *CN* and *<sup>V</sup><sup>m</sup>* stand for the volume fraction of CNTs and matrix, respectively, such that *Vx*<sup>3</sup> *CN* + *Vm* = 1. The density can be defined as

$$V\_{\rm CN}^{\*} = \frac{w\_{\rm CN}}{w\_{\rm CN} + (\rho^{\rm CN}/\rho^{\rm m})(1 - w\_{\rm CN})} \tag{3}$$

whereas the volume fraction for shallow shells takes the following form (see Figure 2)

$$V\_{CN}^{\overline{\pi}\_3} = \begin{cases} \begin{array}{c} UD \ \text{at} \end{array} \begin{array}{c} V\_{CN}^\* \\\\ V D \ \text{at} \ (1-\overline{\pi}\_3) V\_{CN}^\* \\\\ OD \ \text{at} \ (1+\overline{\pi}\_3) V\_{CN}^\* \\\\ \text{XD} \ \text{at} \ \text{4} \ |\overline{\pi}\_3| \ V\_{CN}^\* \end{array} \end{cases} \tag{4}$$

The strain field on the reference surface is governed by the following kinematic relations [36]

$$\begin{array}{ll} \varepsilon\_{11} = \frac{\partial u}{\partial \mathbf{x}\_{1}} - \frac{\overline{w}}{\overline{\mathbf{R}}\_{1}} + \frac{1}{2} \left(\frac{\partial w}{\partial \mathbf{x}\_{1}}\right)^{2}, & \varepsilon\_{22} = \frac{\partial v}{\partial \mathbf{x}\_{2}} - \frac{\overline{w}}{\overline{\mathbf{R}}\_{2}} + \frac{1}{2} \left(\frac{\partial w}{\partial \mathbf{x}\_{2}}\right)^{2} \\\\ \gamma\_{12} = \frac{\partial v}{\partial \mathbf{x}\_{1}} + \frac{\partial u}{\partial \mathbf{x}\_{2}} + \frac{\partial w}{\partial \mathbf{x}\_{1}} \frac{\partial w}{\partial \mathbf{x}\_{2}} \end{array} \tag{5}$$

and the constitutive relations accounting for the von Karman nonlinearity within a classical shell framework are defined as [21]

$$
\begin{bmatrix} \tau\_{11} \\\\ \tau\_{22} \\\\ \tau\_{12} \end{bmatrix} = \begin{bmatrix} E\_{11}^{\overline{\Upsilon}\_{3}} & E\_{12}^{\overline{\Upsilon}\_{3}} & 0 \\\\ E\_{21}^{\overline{\Upsilon}\_{3}} & E\_{22}^{\overline{\Upsilon}\_{3}} & 0 \\\\ 0 & 0 & E\_{66}^{\overline{\Upsilon}\_{3}} \end{bmatrix} \begin{bmatrix} e\_{11} - \overline{\pi}\_{3} \frac{\partial^{2} w}{\partial x\_{1}^{2}} \\\\ e\_{22} - \overline{\pi}\_{3} \frac{\partial^{2} w}{\partial x\_{2}^{2}} \\\\ \gamma\_{12} - 2\overline{\pi}\_{3} \frac{\partial^{2} w}{\partial x\_{1} \partial x\_{2}} \end{bmatrix} \tag{6}
$$

with

$$E\_{\rm ii}^{\mathbb{Z}\_3} = \frac{\mathcal{Y}\_{\rm ii}^{\mathbb{Z}\_3}}{1 - \nu\_{\rm ij}\nu\_{\rm ji}}, \quad E\_{\rm ij}^{\mathbb{Z}\_3} = \frac{\nu\_{\rm ji}\mathcal{Y}\_{\rm ii}^{\mathbb{Z}\_3}}{1 - \nu\_{\rm ij}\nu\_{\rm ji}} = E\_{\rm ji}^{\mathbb{Z}\_3}, \; E\_{66}^{\mathbb{Z}\_3} = \mathcal{Y}\_{\rm ij}^{\mathbb{Z}\_3}(\mathbf{i}, j = 1, 2) \tag{7}$$

**Figure 1.** (**a**) Spherical and (**b**) hypar shallow shells reinforced with CNTs, resting on elastic foundations.

**Figure 2.** Cross-section of shallow shells with uniform and linearly patterned CNTs (**a**) UD, (**b**) VD, (**c**) OD and (**d**) XD.

#### *2.2. Nonlinear Structural Model in the Presence of a PF*

By using relations (1), (2), (5) and (6), the nonlinear governing equations for doubly curved shallow shells reinforced with a linear pattern of CNTs and resting on a WPF, attain the following form

$$L\_{11}(F) + L\_{12}(w) + L\_{13}(F, w) + K(w) = 0\tag{8}$$

$$L\_{21}(F) + L\_{22}(w) + L\_{13}(w, w) = 0\tag{9}$$

where *F* is a stress function and *Lij*(*i* = 1, 2, *j* = 1, 2, 3) are differential operators defined as

*L*11(*F*) = *h u*<sup>12</sup> *∂*4 *∂x*<sup>4</sup> 1 <sup>+</sup> (*u*<sup>11</sup> <sup>−</sup> <sup>2</sup>*u*<sup>31</sup> <sup>+</sup> *<sup>u</sup>*22) *<sup>∂</sup>*<sup>4</sup> *∂x*<sup>2</sup> 1*∂x*<sup>2</sup> 2 + *u*<sup>21</sup> *∂*4 *∂x*<sup>4</sup> 2 + - 1 *R*2 *∂*2 *∂x*<sup>2</sup> 1 + <sup>1</sup> *R*1 *∂*2 *∂x*<sup>2</sup> 2 , *L*12(*w*) = −*u*<sup>13</sup> *∂*4 *∂x*<sup>4</sup> 1 <sup>−</sup> (*u*<sup>14</sup> <sup>+</sup> <sup>2</sup>*u*<sup>32</sup> <sup>+</sup> *<sup>u</sup>*23) *<sup>∂</sup>*<sup>4</sup> *∂x*<sup>2</sup> 1*∂x*<sup>2</sup> 2 − *u*<sup>24</sup> *∂*4 *∂x*<sup>4</sup> 1 − *ρ*<sup>1</sup> *∂*2 *∂t*2 , *L*13(*F*, *w*) = *h* - *∂*2 *∂x*<sup>2</sup> 2 *∂*2 *∂x*<sup>2</sup> 1 <sup>−</sup> <sup>2</sup> *<sup>∂</sup>*<sup>2</sup> *∂x*1*∂x*<sup>2</sup> *∂*2 *<sup>∂</sup>x*1*∂x*<sup>2</sup> <sup>+</sup> *<sup>∂</sup>*<sup>2</sup> *∂x*<sup>2</sup> 1 *∂*2 *∂x*<sup>2</sup> 2 , *K*(*w*) = *kw* − *kp* - *∂*2 *∂x*<sup>2</sup> 1 + *<sup>∂</sup>*<sup>2</sup> *∂x*<sup>2</sup> 2 *L*21(*F*) = *h q*<sup>11</sup> *∂*4 *∂x*<sup>4</sup> 2 + (*q*<sup>12</sup> <sup>+</sup> *<sup>q</sup>*<sup>21</sup> <sup>+</sup> *<sup>q</sup>*31) *<sup>∂</sup>*<sup>4</sup> *∂x*<sup>2</sup> 1*∂x*<sup>2</sup> 2 + *q*<sup>22</sup> *∂*4 *∂x*<sup>4</sup> 1 , *L*22(*w*) = −*q*<sup>23</sup> *∂*4 *∂x*<sup>4</sup> 1 <sup>−</sup> (*q*<sup>24</sup> <sup>+</sup> *<sup>q</sup>*<sup>13</sup> <sup>−</sup> *<sup>q</sup>*32) *<sup>∂</sup>*<sup>4</sup> *∂x*<sup>2</sup> 1*∂x*<sup>2</sup> 2 − *q*<sup>14</sup> *∂*4 *∂x*<sup>4</sup> 2 + - 1 *R*2 *∂*2 *∂x*<sup>2</sup> 1 + <sup>1</sup> *R*1 *∂*2 *∂x*<sup>2</sup> 2 , *L*23(*w*, *w*) = − *∂*<sup>2</sup> *∂x*1*∂x*<sup>2</sup> 2 + *<sup>∂</sup>*<sup>2</sup> *∂x*<sup>2</sup> 1 *∂*2 *∂x*<sup>2</sup> 2 (10)

being *<sup>t</sup>* the time variable, and *<sup>ρ</sup>*<sup>1</sup> <sup>=</sup> *<sup>h</sup>* /2 −*h*/2 *ρx*3 <sup>1</sup> *dx*3. Moreover, *uij* are defined as

$$\begin{aligned} u\_{11} &= A\_{11}^1 q\_{11} + A\_{12}^1 q\_{21}, \; u\_{12} = A\_{11}^1 q\_{12} + A\_{12}^1 q\_{11}, \; u\_{13} = A\_{11}^1 q\_{13} + A\_{12}^1 q\_{23} + A\_{11}^2, \\ u\_{14} &= A\_{11}^1 q\_{14} + A\_{12}^1 q\_{24} + A\_{12}^2, \; u\_{21} = A\_{21}^1 q\_{11} + A\_{22}^1 q\_{21}, \; u\_{22} = A\_{21}^1 q\_{12} + A\_{22}^1 q\_{22}, \\ u\_{23} &= A\_{21}^1 q\_{13} + A\_{22}^1 q\_{23} + A\_{21}^1, \; u\_{24} = A\_{21}^1 q\_{14} + A\_{22}^1 q\_{24} + A\_{22}^1, \; u\_{31} = A\_{66}^1 q\_{35}, \\ u\_{32} &= A\_{66}^1 q\_{32} + 2 A\_{66}^2. \end{aligned} \tag{11}$$

with

$$\begin{aligned} q\_{11} &= \frac{A\_{22}^0}{\Pi'}, q\_{12} = -\frac{A\_{12}^0}{\Pi'}, q\_{13} = \frac{A\_{12}^0 A\_{21}^1 - A\_{11}^0 A\_{22}^0}{\Pi}, q\_{14} = \frac{A\_{12}^0 A\_{21}^2 - A\_{12}^0 A\_{22}^0}{\Pi}, q\_{21} = -\frac{A\_{21}^0}{\Pi'}, q\_{22} = \frac{A\_{11}^0}{\Pi'},\\ q\_{23} &= \frac{A\_{11}^1 A\_{21}^0 - A\_{21}^1 A\_{11}^0}{\Pi}, q\_{24} = \frac{A\_{12}^1 A\_{21}^0 - A\_{22}^1 A\_{11}^0}{\Pi}, q\_{31} = \frac{1}{A\_{66}^0}, q\_{32} = -\frac{2A\_{66}^1}{A\_{66}^0} \Pi = A\_{11}^0 A\_{22}^0 - A\_{12}^0 A\_{21}^0,\\ A\_{ij}^{i\_1} &= \int\_{\bf{j}}^{h/2} \underline{E}\_{k\bf{k}}^{\sf T} \underline{\bf{j}}^{\bf j} \mathrm{d} \overline{\bf{x}}\_{3} (k = 1, 2, 6; \, i\_1 = 0, 1, 2), \quad A\_{ij}^{i\_1} = \int\_{-h/2}^{h/2} \underline{\bf{x}}\_{i\bf{j}}^{\sf T} \underline{\bf{x}}\_{3}^{\bf j} \mathrm{d} \overline{\bf{x}}\_{3} (i, j = 1, 2). \end{aligned} \tag{12}$$

#### **3. Solution Procedure**

In what follows, we provide an analytical solution to the problem of a simplysupported doubly-curved shell. Thus, the structural deflection can be approximated as [21,36]

$$w = \overline{w}(t)\sin(\mathfrak{a}\_1\mathfrak{x}\_1)\sin(\mathfrak{a}\_2\mathfrak{x}\_2) \tag{13}$$

where *w*(*t*) is a function of time, *α*<sup>1</sup> = *<sup>m</sup><sup>π</sup> <sup>a</sup>* , *<sup>α</sup>*<sup>2</sup> <sup>=</sup> *<sup>n</sup><sup>π</sup> <sup>b</sup>* , in which *m* and *n* are the wave numbers in directions *x*<sup>1</sup> and *x*2, respectively. By substitution of Equation (13) into Equation (9), we get the following expression for the stress function *F*

$$F = \overline{w}(t)[c\_1\overline{w}(t)\cos(2u\_1\mathbf{x}\_1) + c\_2\overline{w}(t)\cos(2u\_2\mathbf{x}\_2) + c\_3\sin(u\_1\mathbf{x}\_1)\sin(u\_2\mathbf{x}\_2)] \tag{14}$$

where *cj*(*j* = 1, 2, 3) are defined as

$$c\_1 = \frac{a\_2^2}{32a\_1^2 q\_2 n'}, \ c\_2 = \frac{a\_1^2}{32a\_2^2 q\_{11} h'}, \ c\_3 = \frac{q\_{23} a\_1^4 + (q\_{24} + q\_{13} - q\_{32})a\_1^2 a\_2^2 + q\_{14} a\_2^4 + a\_1^2 / R\_2 + a\_2^2 / R\_1}{h \left[ q\_{11} a\_2^3 + (q\_{12} + q\_{21} + q\_{31})a\_1^2 a\_2^2 + q\_{22} a\_1^4 \right]} \tag{15}$$

By substituting Equations (13) and (15) into Equation (8) and by applying the Galerkin procedure in the domain 0 ≤ *x*<sup>1</sup> ≤ *a* and 0 ≤ *x*<sup>2</sup> ≤ *b*, we obtain

$$L(t) \equiv \frac{d^2 \tilde{w}(t)}{dt^2} + \left(\omega\_{wp}^L\right) \tilde{w}(t) + \theta\_1 \tilde{w}^2(t) + \theta\_2 \tilde{w}^3(t) = 0 \tag{16}$$

where *<sup>w</sup>* (*t*) = *<sup>w</sup>*(*t*)/*h*, the quantities *<sup>θ</sup>*1, *<sup>θ</sup>*<sup>2</sup> are defined as

$$\theta\_1 = \frac{64b^2}{34b} \frac{1}{\rho\_1} \left( \frac{\mu\_1 a\_1^4 \epsilon\_1 + \mu\_2 a\_2^4 \epsilon\_2}{a\_2 a\_1} - \frac{a\_1 a\_2^\* \epsilon\_2}{8} - \frac{c\_1}{4\overline{k}\_2} \frac{a\_1}{a\_2} - \frac{c\_2}{4\overline{k}\_1} \frac{a\_2}{a\_1} \right) \left[ 1 - (-1)^m - (-1)^n + (-1)^{m+n} \right] \tag{17}$$
 
$$\theta\_2 = \frac{2b^3 a\_1^2 a\_2^2 (\epsilon\_1 + \epsilon\_2)}{\rho\_1}$$

and *ω<sup>L</sup> wp* is the frequency associated to the shallow structure resting on the PF at small deflections, defined as

$$\begin{aligned} \omega\_{wp}^L &= \frac{1}{\sqrt{\rho\_1}} \left\{ \left[ \frac{a\_1^2}{R\_2} + \frac{a\_2^2}{R\_1} - \mu\_{12} a\_1^4 - (\mu\_{11} - 2\mu\_{31} + \mu\_{22}) a\_1^2 a\_2^2 - \mu\_{21} a\_2^4 \right] \|s\_3\\ &+ \mu\_{13} a\_1^4 + (\mu\_{14} + 2\mu\_{32} + \mu\_{24}) a\_1^2 a\_2^2 + \mu\_{24} a\_2^4 + k\_w + k\_p (a\_1^2 + a\_2^2) \right\}^{1/2} \end{aligned} \tag{18}$$

The approximate solution of Equation (16) reads as follows

$$
\tilde{w}(t) = w\_0 \cos(\mathcal{O}\_{NL} t) \tag{19}
$$

where *w*<sup>0</sup> is the dimensionless amplitude, *NL* is the nonlinear frequency and the initial conditions are defined as

$$\left. \tilde{w}(t) \right|\_{t=0} = w\_0 \quad \text{and} \left. \frac{\text{d}\tilde{w}(t)}{\text{d}t} \right|\_{t=0} = 0 \tag{20}$$

By combining the relations (16) and (19), we obtain an equation of the type *L*(*t*) = 0. Thus, by applying the Grigolyuk method [37], one obtains

$$\int\_{0}^{\pi/2\mathfrak{a}\_{NL}} L(t) \cos(\mathfrak{o}\_{NL}t)dt = 0\tag{21}$$

After integrating this last relation, we obtain the following nonlinear amplitude– frequency dependence

$$
\omega\_{wp}^{NL} = \sqrt{\left(\omega\_{wp}^{L}\right)^2 + \frac{8}{3\pi}\theta\_1 w\_0 + 0.75\theta\_2 w\_0^2} \tag{22}
$$

and

$$
\boldsymbol{\phi}\_{1wp}^{NL} = \boldsymbol{\phi}\_{wp}^{NL} \boldsymbol{h} \sqrt{\boldsymbol{\rho}^m/\boldsymbol{Y}^m} \tag{23}
$$

The rational nonlinear-to-linear free vibration frequency (NLFVF / LFVF), *NL wp* /*ωL*, becomes

$$\frac{\omega\_{wp}^{NL}}{\omega\_L} = \sqrt{\left(\frac{\omega\_{wp}^L}{\omega\_L}\right)^2 + \frac{8}{3\pi} \frac{\theta\_1 w\_0}{\omega\_L^2} + \frac{0.75 \theta\_2 w\_0^2}{\omega\_L^2}}\tag{24}$$

where the linear-free vibration frequency *ω<sup>L</sup>* for the unconstrained structure is defined as (18) and *kw* = *kp* = 0 represents a special case. Based on Equations (22) and (24), we can treat different cases, namely shallow spherical shells (for *R*<sup>1</sup> = *R*2) or shallow hypar shells (for *R*<sup>1</sup> = −*R*2) resting on a PF, as well as shallow cylindrical panels (*R*<sup>1</sup> → ∞) or plates (*R*<sup>1</sup> → ∞, *R*<sup>2</sup> → ∞), on a WPF.

The lowest values of *NL wp* , *NL* <sup>1</sup>*wp* and *NL wp <sup>ω</sup><sup>L</sup>* for shallow spherical and hypar shells on a PF are determined by minimizing Equations (20)–(22) depending on the vibration modes (*m, n*), for fixed values of the dimensionless amplitude *w*0.

#### **4. Results and Discussion**

The numerical investigation starts with a comparative evaluation of the dimensionless linear frequency parameters, *ω*1*<sup>L</sup>* = *ωLh* /*ρm*/*Ym*, for isotropic shallow shells with respect to predictions from the literature [38] (see Table 1).

**Table 1.** Comparative evaluation with the literature of *ω*1*<sup>L</sup>* for different shallow structural members made of an isotropic material.


Thus, we use the expression (18), while keeping *kw* = *kp* = 0, *Y<sup>m</sup>* <sup>11</sup> = *<sup>Y</sup><sup>m</sup>* <sup>22</sup> = *<sup>Y</sup><sup>m</sup>* = 70 Gpa, *<sup>ν</sup>*<sup>12</sup> = *<sup>ν</sup>*<sup>21</sup> = *<sup>ν</sup><sup>m</sup>* = 0.3177,*ρ<sup>m</sup>* = 2.702 × 103 kg/m<sup>3</sup> and the geometrical ratios a/b = 1, a/h = 10. As visible from Table 1, our results match very well predictions from [38], for both spherical and hypar shell members; this proves the reliability and consistency of the proposed formulation. A further comparison with the literature [39,40] is also provided in terms of dimensionless linear frequencies for isotropic square plates resting on PFs with *h*/*b* = 0.01, *a*/*b* = 1, *kw* = 100*D<sup>m</sup>* and *kp* = 10*Dm*. For this subcase, the relation (16) is computed for *R*<sup>1</sup> → ∞, *R*<sup>2</sup> → ∞ and the dimensionless linear frequency parameter is determined as Ω*Lwp* = *ω<sup>L</sup> wp*(*b*/*π*) <sup>2</sup>/*ρmh*/*D<sup>m</sup>* in which *D<sup>m</sup>* = *<sup>Y</sup>mh*<sup>3</sup> <sup>12</sup>[1−(*νm*) 2 ] , see [40]. Table 2 summarizes the results based on different approaches, with a consistent agreement between our formulation and findings from [39,40].

**Table 2.** Comparison of dimensionless frequency parameters for square plates on a PF (*h*/*b* = 0.01, *a*/*b* = 1, *kw* = 100*Dm*, *kp* = 10*D<sup>m</sup>* ).


After this preliminary validation, we continue the analysis by computing the NLFVF of shallow spherical and hypar shells reinforced with a uniform and linear distribution of CNTs, and resting on a WF and PF. The selected shell members are made of polymethyl methacrylate (PMMA), as matrix, and single-walled CNTs, with geometrical properties *<sup>r</sup>*<sup>1</sup> = 9.26 nm, *<sup>a</sup>*<sup>1</sup> = 6.8 × <sup>10</sup>−<sup>1</sup> nm, *<sup>h</sup>*<sup>1</sup> = 6.7 × <sup>10</sup>−<sup>2</sup> nm, as reinforcement. The mechanical properties for the CNT phase are *YCN* <sup>11</sup> = 5646.6 Gpa, *<sup>Y</sup>CN* <sup>22</sup> = <sup>7080</sup> Gpa, *<sup>Y</sup>CN* <sup>12</sup> = 1944.5 Gpa, *νCN* <sup>12</sup> = 0.175, *<sup>ρ</sup>CN* = 1.4 × 103 kg/m3; for PMMA, it is *<sup>Y</sup><sup>m</sup>* = 2.5 Gpa, *<sup>ν</sup><sup>m</sup>* = 0.34, *<sup>ρ</sup><sup>m</sup>* = 1.15 × <sup>10</sup><sup>3</sup> kg/m3. In line with [9], we also consider different efficiency parameters of CNT/matrix depending on the selected value of *V*∗ *CN*, as summarized in Table 3. As also listed in Table 4, we check for the variation of NLFVFs for both the selected shallow shells, while keeping different distributions of CNTs (i.e., UD, VD, OD and XD), and by varying the stiffness constants (*kw*, *kp*) for the elastic foundation under the three fixed values of *V*∗ *CN* (0.12, 0.17, and 0.28). The frequency values are computed for *R*<sup>1</sup> = 20*h*, *a* = *b*, *a* = 20*h*, (*m*, *n*)=(1, 1) and *w*<sup>0</sup> = 1.5, with a clear increase in results for an increased value of the stiffness parameters, *kp* and/or *kw*, for all the reinforcement assumptions. For fixed

values of *kp*, *kw*, and *V*∗ *CN*, it also seems that OD and XD patterns of CNTs always provide the lowest and highest frequency values, respectively, independently of the selected shell geometry.

**Table 3.** Typical properties of CNT/matrix efficiency parameters depending on the volume fraction of CNTs.


**Table 4.** Variation in NLFVF for shallow spherical and hypar shells with CNTs resting on a W-EF and PF with various foundation elastic parameters, versus *V*∗ *CN*.


The highest sensitivity of the response to the CNT dispersion within the matrix is observed for a fixed value of *V*∗ *CN* = 0.28, with a maximum percentage variation with respect to a UD of 21.6%. At the same time, the largest foundation effect on *NL* <sup>1</sup>*wp* occurs at *V*∗ *CN* = 0.12 and OD patterns with a percentage variation of 51.5%. A PF also seems to affect the response more significantly compared to a W-EF, reaching the highest sensitivity with an OD-type reinforcement and *V*∗ *CN* = 0.12, whereas the lowest sensitivity is obtained for a XD of CNTs with *V*∗ *CN* = 0.28.

As far as the sensitivity to the volume fraction is concerned, the largest influence is noticed for structures on a PF reinforced by XD CNTs with *V*∗ *CN* = 0.28 and the lowest

effect is obtained with an OD pattern and *V*∗ *CN* = 0.12, respectively. In Figure 3, we plot the variation in NLFVFs for shallow spherical and hypar shells reinforced with a UD and VD of CNTs versus *w*0, for a fixed value of *V*∗ *CN* = 0.28, for three different geometrical ratios *R*1/*a* = 2.0, 2.5, 3.0, accounting (or not) for the presence of a surrounding WPF. The other parametric data are: *<sup>a</sup>*/*<sup>b</sup>* = 0.5, *<sup>a</sup>*/*<sup>h</sup>* = 15, (*m*, *<sup>n</sup>*)=(1, 1), *kw* = <sup>4</sup> × <sup>10</sup>9(N/m3) and *kp* = 1.6 × <sup>10</sup>4(N/m). As visible in Figure 3, the NLFVF for hypar shells resting on a PF increases monotonically with *w*0, whereas it varies non-monotonically for spherical shells on a PF with an initial decrease for *w*<sup>0</sup> ≤ 0.5, and a further increase for *w*<sup>0</sup> > 0.5. By comparing results among unconstrained spherical and hypar members, it seems that NLFVFs for hypar shells always reach higher values than spherical ones for all *w*0; the NLFVFs for shallow spherical shells are usually higher for *w*<sup>0</sup> ≤ 0.5.

**Figure 3.** Variation in NLFVF for shallow spherical and hypar shells with UD- and VD-patterned CNTs on a PF versus *w*0, with different geometrical ratios *R*1/*a*.

The magnitude of NLFVFs for shell members with and without a PF can vary significantly under the same geometrical assumption *R*1/*a* and the same value of *w*0. In addition, for an increasing rational value of *R*1/*a*, the NLFVF values of spherical shells decrease for *w*<sup>0</sup> ≤ 0.5 and increase for *w*<sup>0</sup> > 0.5. The influence of CNT patterns on the NLFVF of hypar shells in the presence, or not, of a PF, decreases with a varying *w*0. For spherical shells, instead, such an effect increases for *w*<sup>0</sup> ≤ 0.5 and decreases for *w*<sup>0</sup> > 0.5. Such sensitivity becomes more pronounced for *w*<sup>0</sup> > 0.5. More specifically, the influence of CNT patterns on the NLFVF for unconstrained hypar shells decreases from −16.71% to −5.41% due to the increase in *w*0. For unconstrained spherical shells, the influence of CNT patterns increases from −15.20% to −17.51% in the range of *w*<sup>0</sup> ≤ 0.5, and decreases from −17.51% to −10.82% for *w*<sup>0</sup> > 0.5, under a fixed ratio *R*1/*a* = 2. The influence of different CNT patterns is more pronounced for unconstrained hypar shells, with the largest difference being approximately 1.00%; for spherical shells on a PF, the largest difference becomes approximately 1.8% due to an increased ratio *R*1/*a*.

The influence of CNT patterns on NLFVFs reduces with a maximum percentage of 4.76% and 4.61%, for spherical and hypar shells on PF, respectively. A pronounced effect of CNT patterns is also observed for both shallow shells in the presence, or not, of a PF, which is quantified as a percentage by 4.63% and 7.19%, respectively. The effect of a PF on NLFVF for both shells is approximately 6% greater for a VD pattern compared to a UD pattern.

Variation in the NLFVF / LFVF ratio, for shallow spherical and hypar shells with UD and OD patterns versus *w*0, is plotted in Figure 4, for *V*∗ *CN* = 0.12; 0.17; 0.28, while keeping *<sup>a</sup>*/*<sup>h</sup>* = <sup>10</sup> *<sup>R</sup>*1/*<sup>a</sup>* = 2, *<sup>a</sup>*/*<sup>b</sup>* = 2, (*m*, *<sup>n</sup>*)=(1, 1), *kw* = <sup>3</sup> × <sup>10</sup>9(N/m3) and *kp* = 1.5 × 104(N/m). As clearly visible in Figure 4, the NLFVF / LFVF ratio varies nonmonotonically for shallow spherical shells and monotonically for hypar shells, with an increased value of *w*<sup>0</sup> for all selected values of *V*∗ *CN*.

**Figure 4.** Variation in the NLFVF / LFVF ratio for shallow spherical and hypar shells with UD- and OD-patterned CNTs, in the presence/absence of a PF versus *w*0, with different *V*∗ *CN*.

Based on a comparison of the NLFVF / LFVF ratio for hypar shells with UD and OD patterns, in presence or absence of a PF, a higher variation is noticed for an OD pattern. More specifically, the NLFVF / LFVF ratio for spherical shells with an OD pattern becomes higher for all *w*<sup>0</sup> in presence of a PF, and for *w*<sup>0</sup> > 0.5 in the absence of a surrounding elastic medium. When the NLFVF / LFVF ratio is evaluated comparatively for spherical shells with different *V*∗ *CN*, the largest NLFVF / LFVF ratio for unconstrained spherical shells occurs at *V*∗ *CN* = 0.28, and for spherical shells on PF at *V*<sup>∗</sup> *CN* = 0.12. For hypar shells in presence or not of a surrounding elastic medium, the highest NLFVF / LFVF ratio is always obtained for *V*∗ *CN* = 0.28. It is also noticeable that this ratio becomes higher for hypar shells with and without the PF, as spherical and hypar shells are compared. The pattern effect on the NLFVF / LFVF ratio increases for unconstrained spherical shells (4%) and unconstrained hypar shells (5%), whereas the pattern effect on the *NL wp* /*ω<sup>L</sup>* ratio in both shells on PF increases with *w*0, accordingly. The most pronounced increase seems to be approximately equal to 2.5% for spherical shells, and approximately equal to 2.1% for hypar shells. In absence of a surrounding elastic medium, the influence of a CNT pattern on *NL wp* /*ω<sup>L</sup>* for hypar shells becomes 1.9% more pronounced than unconstrained spherical shells on a PF.

The effect of a PF on *NL wp* /*ω<sup>L</sup>* ratio decreases for an increased value of *w*0, for both spherical and hypar shells, with a maximum percentage variation of 6% and 10%, respectively.

The variation in *NL wp* /*ω<sup>L</sup>* with the PF is about 3.5% for spherical shells with an OD pattern of CNTs; this effect is 2.4% more pronounced for hypar shells. In addition, for a reinforcement phase with *V*∗ *CN* = 0.12, the percentage variation of *NL wp* /*ω<sup>L</sup>* for both shells under a PF is approximately 6% (or 9%) greater than that one for *V*∗ *CN* = 0.17 (or *V*∗ *CN* = 0.28).

Figure 5 shows the variation in the NLFVF / LFVF ratio of spherical shells on PF, reinforced with UD- and XD-patterned CNTs, against *w*0, for *V*∗ *CN* = 0.17. In this parametric study we also consider different values of *a*/*b* (i.e., *a*/*b* = 0.5, 1.0, 1.5), along with *<sup>R</sup>*1/*<sup>a</sup>* = 2, *<sup>a</sup>*/*<sup>h</sup>* = 15, (*m*, *<sup>n</sup>*)=(1, 1), *kw* = <sup>3</sup> × <sup>10</sup>9(Pa/m) and *kp* = 1.5 × 105(Pa.m). As visible in Figure 5, for an increased value of *a*/*b*, the NLFVF / LFVF ratios of spherical shells with and without a PF vary nonmonotonically, with an increase after an initial decrease up to a minimum value. The NLFVF / LFVF ratio for UD patterns is larger than XD patterns for all *w*<sup>0</sup> (in presence of PF) and for *w*<sup>0</sup> > 1 (in absence of an elastic ground). The same ratio, for XD patterns, is larger than UD patterns, as *w*<sup>0</sup> ≤ 1 only in the absence of ground.

**Figure 5.** Variation in the NLFVF / LFVF ratio for shallow spherical shells containing UD- and XD-patterned CNTs in the presence/absence of a PF versus *w*0, with different values of *a*/*b*.

For an increased value of *a*/*b*, the NLFVF / LFVF ratio of unconstrained spherical shells decreases for both patterns, while it decreases (or increases) when *w*<sup>0</sup> > 1.25 (or *w*<sup>0</sup> ≤ 1.25), for spherical shells on a PF. It is also observed that XD patterns effect on the NLFVF / LFVF ratio is higher in presence of a PF; it decreases/increases depending on the value of *w*0, while it continuously decreases depending on the increase in *a*/*b*. Although the effect of PF on NLFVF / LFVF ratio is more pronounced for UD patterns, it decreases depending on the increase in *w*<sup>0</sup> and increases for an increased value of *a*/*b*. The minimum and maximum influence of PF on the NLFVF / LFVF ratio for spherical shells corresponds to a percentage variation of 10.24% and 24.71%, respectively.

Finally, in Figure 6 we plot the variation of NLFVF / LFVF ratio for UD- and VDpatterned spherical and hypar shells (in the presence or absence of a PF), versus *w*<sup>0</sup> for different *R*1/*a* ratios (i.e., *R*1/*a* = 2.0, 2.5, 3.0), while keeping *V*<sup>∗</sup> *CN* = 0.28, *a*/*b* = 0.5, *a*/*h* = 15, (*m*, *<sup>n</sup>*)=(1, 1), *kw* = <sup>3</sup> × <sup>10</sup>9(N/m3) and *kp* = 1.5 × <sup>10</sup><sup>5</sup> (N/m). The NLFVF / LFVF ratio of hypar shells on PF increases with *w*0, while it varies nonmonotonically with *R*1/*a*. Similarly, the NLFVF / LFVF ratio of spherical shells on PF decreases first and then increases for an increased value of *w*0, while always increasing for an increased *R*1/*a* ratio. The NLFVF / LFVF ratio of hypar shells for VD patterns with and without a PF, as well as for spherical shells resting on a PF, is greater than the same shells reinforced uniformly by CNTs. The NLFVF / LFVF ratio of unconstrained spherical shells with VD patterns is higher than those with UD patterns, at least when *w*<sup>0</sup> ≤ 0.75; the contrary occurs for *w*<sup>0</sup> > 0.75. Based on a comparative evaluation of both geometries, the NLFVF / LFVF ratios of UD- and VD-reinforced spherical shells on a PF are lower than hypar shells. The influence of VD patterns on *NL wp* /*ω<sup>L</sup>* for spherical shells is higher than hypar shells, with a maximum increase of 3.2% or 0.6%, respectively, depending on the increase in *R*1/*a*. Looking at the influence of the foundation on *NL wp* /*ω<sup>L</sup>* for spherical shells with a UD of CNTs, it decreases first, up to a minimum value, then increases for an increased value of *w*0. A monotonic decrease is differently observed for hypar shells. Depending on the increase of *R*1/*a*, the effect of PF is lower than 2% for spherical shells, and lower than 1.5% for hypar shells.

**Figure 6.** Variation in the NLFVF/LFVF ratio for shallow spherical and hypar shells with UD- and VD-patterned CNTs, in the presence/absence of a PF versus *w*<sup>0</sup> with different values of *R*1/*a*.

#### **5. Conclusions**

In this work, the Donnell's nonlinear shell theory is applied to study the free vibration behavior of composite shell structures reinforced by uniform and linearly patterned CNTs resting on a PF. Once the basic relations for composite shallow shells reinforced by CNTs on WPFs are established, the partial differential equations of nonlinear motion are derived, taking into account the von Karman nonlinearity. These equations are solved here by means of the Galerkin and Grigolyuk methods in terms of linear and nonlinear free vibrations for inhomogeneous nanocomposite construction members such as plates, panels, spherical and

hyperbolic paraboloidal (hypar) shallow shells. The accuracy of the results in the current study has been confirmed by means of a successful comparison with reliable predictions from the literature. After this preliminary validation, a detailed numerical analysis is performed, including the effect of nonlinearity, CNT patterns and volume fraction on the nonlinear frequency response. Based on a large systematic investigation, the analytical results could serve as valid benchmark solutions for further computational studies on the topic, as well as for design purposes. Among the most useful insights, it is found that the variation rate of NLFVFs for both shallow shells with linearly patterned CNTs decreases, while remaining constant for different elastic foundations with an increased stiffness. For both shallow shells, a single- or dual-parameter elastic foundation yields an increase in NLFVFs, where the NLFVF decreases for VD, OD and XD patterns, as foundation coefficients increase. Moreover, the influence of PF on the *NL* <sup>1</sup>*wp* for shallow spherical and hypar shells reinforced with CNTs has revealed as more pronounced than that of a WEF. The highest influence of PF on NLFVF values is observed with an OD pattern of CNTs for *V*∗ *CN* = 0.12, whereas the smallest effect is observed with an XD pattern of CNTs for *V*∗ *CN* = 0.28, respectively, when the influences of PF on NFVFs for spherical or hypar shells are compared to each other. Based on a comparative evaluation of the nonlinear vibration response for both shallow shells, the largest effect of the PF on NLFVFs is observed for a XD CNT-based reinforcement with a volume fraction, *V*∗ *CN* = 0.28, whereas the smallest effect occurs for an OD pattern and *V*∗ *CN* = 0.12. At the same time, the NLFVF of hypar shells on PF increases continuously for an increased *w*0; it decreases when *w*<sup>0</sup> ≤ 0.5 for spherical shells on PF and increases for *w*<sup>0</sup> > 0.5. The pattern effect on *NL wp* /*ω<sup>L</sup>* ratio for spherical shells (4%) and hypar shells (5%) increases with *w*<sup>0</sup> in absence of a PF. Based on the parametric study, the influence of the PF on *NL wp* /*ω<sup>L</sup>* ratio seems to be more pronounced for a UD pattern, but it decreases depending on the increase in *w*<sup>0</sup> and increases for an increased geometrical ratio *a*/*b*. Moreover, the rational value of *NL wp* /*ω<sup>L</sup>* for a UD pattern is higher than a XD pattern, for all *w*<sup>0</sup> in presence of an elastic medium, and for *w*<sup>0</sup> > 1 in absence of an elastic medium. The same ratio for an XD pattern is larger than the one for a UD pattern, as *w*<sup>0</sup> ≤ 1 only in absence of an elastic ground. Finally, for a VD pattern, the rational value of *NL wp* /*ω<sup>L</sup>* for spherical shells is greater than the hypar shells, with a maximum increase of 3.2% in lieu of 0.6%, as found for hypar shells, depending on the increase in the geometrical ratio *R*1/*a*.

**Author Contributions:** Conceptualization, A.M., F.T., R.D. and N.K.; Formal analysis, A.M., F.T., R.D. and N.K.; Investigation, A.M. and N.K.; Validation, A.M., F.T., R.D. and N.K.; Writing—Original Draft, A.M. and N.K.; Writing—Review & Editing, F.T. and R.D. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research received no external funding.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


## *Article* **Forced Vibration Analysis of Composite Beams Reinforced by Carbon Nanotubes**

**Ömer Civalek 1, ¸Seref D. Akba¸s 2, Bekir Akgöz <sup>3</sup> and Shahriar Dastjerdi 3,4,\***


**Abstract:** This paper presents forced vibration analysis of a simply supported beam made of carbon nanotube-reinforced composite material subjected to a harmonic point load at the midpoint of beam. The composite beam is made of a polymeric matrix and reinforced the single-walled carbon nanotubes with their various distributions. In the beam kinematics, the first-order shear deformation beam theory was used. The governing equations of problem were derived by using the Lagrange procedure. In the solution of the problem, the Ritz method was used, and algebraic polynomials were employed with the trivial functions for the Ritz method. In the solution of the forced vibration problem, the Newmark average acceleration method was applied in the time history. In the numerical examples, the effects of carbon nanotube volume fraction, aspect ratio, and dynamic parameters on the forced vibration response of carbon nanotube-reinforced composite beams are investigated. In addition, some comparison studies were performed, with special results of published papers to validate the using formulations.

**Keywords:** carbon nanotube-reinforced composite; forced vibration; dynamic analysis; beam; harmonic load

#### **1. Introduction**

Composite material refers to any solid that consists of more than one component, in which they are in separate phases. The main advantages of composite materials are excellent strength-to-weight and stiffness-to-weight ratios. The fibrous composites, consisting of carbon, glass, aramid, and basalt fibers, have a wide range of applications in many modern engineering and industries, such as civil, automotive, bicycle, mechanical, defense, marine, aviation, and aerospace [1–9].

In the early 1990s, carbon nanotubes (CNTs) were studied by Sumio Iijima [10]. CNTs can be considered as one of the key building blocks of nanotechnology. Interest of the scientists and researchers in CNTs' potential engineering and its industrial applications, such as in aerospace, composites, electronics, computers, energy, medicine, sensors, and air and water purifications, have grown rapidly due to the unique material properties of CNTs [11–18].

Due to extraordinary electrical, thermal, and mechanical properties, besides providing good interfacial bonds, CNTs can be ideally suited for the next generation of composite materials and have recently been used instead of the traditional fibers for the reinforcing element in the reinforced composites. Some studies on the bending, buckling, and vibration responses of carbon nanotube-reinforced composite (CNTRC) beams are mentioned below.

Free vibration analysis of functionally graded single-walled CNTs reinforced aluminum alloy beam was performed by Selmi and Bisharat [19]. They obtained the natural

**Citation:** Civalek, Ö.; Akba¸s, ¸S.D.; Akgöz, B.; Dastjerdi, S. Forced Vibration Analysis of Composite Beams Reinforced by Carbon Nanotubes. *Nanomaterials* **2021**, *11*, 571. https://doi.org/10.3390/nano 11030571

Academic Editor: Raffaele Barretta

Received: 2 February 2021 Accepted: 22 February 2021 Published: 25 February 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

frequencies by both the Rayleigh–Ritz method and ANSYS, a software package that implements the finite element method. Pure bending and local buckling behaviors of a nanocomposite beam reinforced by a single-walled carbon nanotube were investigated by Vodenitcharova and Zhang [20]. Ke et al. [21] surveyed the nonlinear free vibration response of functionally graded carbon nanotube-reinforced composite (FG-CNTRC) beams on the basis of first-order shear deformation beam theory with von Kármán geometric nonlinearity assumption. The Ritz method with a direct iterative method was used to obtain solutions for different boundary conditions. Yas and Heshmati [22] investigated the vibrational characteristics of FG-CNTRC beams reinforced by randomly oriented CNTs subjected to a moving load, based on Bernoulli–Euler and Timoshenko beam theories by finite element method. Yas and Samadi [23] analyzed the free vibration and buckling problems of embedded nanocomposite beams resting on an elastic foundation. SWCNTs were used as reinforcement elements. Four different CNTs distributions along the height were taken into consideration. The governing equations were derived by implementing Hamilton's principle, and the resulting equations were solved by generalized differential quadrature method. Deepak et al. [24] developed a spectral finite element formulation for uniform and tapered rotating CNTRC polymer beams. The obtained results were comparatively presented with the results of carbon fiber-reinforced laminated composite rotating beams. Ke et al. [25] examined the dynamic stability analysis of FG-CNTRC beams subjected to axial loading. The related governing equations were derived based on Timoshenko beam theory and they were solved by differential quadrature method.

Heshmati and Yas [26] perused the free vibration characteristics of FG-CNTRC beams by using an equivalent fiber based on the Eshelby–Mori–Tanaka approach via the finite element method. The influences of agglomeration and distribution of carbon nanotubes were investigated in detail. Wattanasakulpong and Ungbhakorn [27] performed the static bending, free vibration, and buckling analyses of embedded CNTRC beams lying on the Pasternak elastic foundation. Analytic solutions were achieved by Navier's solution technique for simply supported CNTRC beams. Linear and nonlinear vibrations of CNTRC beams were comprehensively studied on the basis of first and third-order shear deformation beam theories [28–32]. Mayandi and Jeyaraj [33] investigated the static and dynamic behaviors of FG-CNTRC beams subjected to various non-uniform thermal loads by employing the finite element method. Fattahi and Safaei [34] performed stability analysis of CNTRC beams based on three different beam theories. The obtained equations were solved by generalized differential quadrature method. Babu Arumugam et al. [35] surveyed the free and forced vibration responses of CNTRC beams with constant and variable cross-sections by using the finite element method. A detailed parametric work was made to show the influences of slenderness ratio, percentage of CNT constituent, ply orientation, and boundary conditions. Mohseni and Shakouri [36] studied the free vibration and buckling responses of tapered FG-CNTRC beams surrounding an elastic medium based on Timoshenko beam theory. Dynamic analysis of the pre-twisted FG-CNTRC beams subjected to thermal loading was carried out by Shenas et al. [37]. A higher-order shear deformation beam theory was employed to derive the constitutive equations, and the Chebyshev–Ritz method was applied to solve the resulting equations for various end conditions. Khosravi et al. [38] performed the thermal stability analysis of rotating CNTRC beams. Timoshenko beam theory was used in the derivation of the governing differential equations. Generalized differential quadrature method was employed to obtain some numerical results. Additionally, buckling and vibration responses of micro-and nanocomposite beams were investigated [39–41]. On the other hand, there have recently been many studies on the mechanical behaviors of CNTRC plates and shells [42–56].

As mentioned above, such structures may be subjected to various loads. Due to this fact, it is very important to understand the dynamical behavior of these structures subjected to the harmonic loads. Additionally, the aforementioned review reveals that researchers have so far examined the bending, buckling, and free vibration responses of CNTRC beams. In particular, to the best of the authors' knowledge, forced vibration of CNTRC beams due

to harmonic loads has not been investigated in detail. The main purpose of the present study is to fill this gap.

In this paper, forced vibrational behavior of CNTRC beam is examined. It was considered that CNTRC beams are made of a polymeric matrix reinforced by the single-walled carbon nanotubes and is subjected to a harmonic point load at the middle. Three different distributions of CNTs are considered in the analysis. The governing equations of problem have been derived by using the Lagrange procedure based on Timoshenko beam theory. The Ritz method, in conjunction with algebraic polynomials with the trivial functions, was utilized to solve the resulting equation. Additionally, the Newmark average acceleration method was used in the time history for the solution of the forced vibration problem. A detailed parametric study was carried out to peruse the influences of CNTs volume fraction, slenderness ratio, and dynamic parameters on the forced vibration response of CNTRC beam.

#### **2. Theory and Formulation**

Consider a simply supported beam made of CNTRC material under a dynamic load, as shown in Figure 1. The composite beam was subjected to a dynamic point load, Q(t), at the midpoint. The geometry of the beam was indicated as the length, *L*; the height, *h*; and width, *b*. Additionally, three different patterns of CNTs, such as uniformly distributed (UD), Λ- and X- type distributions, were considered throughout the thickness of the composite beam.

**Figure 1.** A schematic illustration of a simply supported carbon nanotube-reinforced composite (CNTRC) beam under a harmonic load and three different patterns of CNTs.

The dynamic point load (*Q*(*t*)) was assumed to be sinusoidal harmonic in time domains, such as in the following:

$$\mathbf{x}(t) = \mathbb{Q}\_0 \sin(\overline{\omega}t), \; 0 \le t \ll \infty \tag{1}$$

where *Q*<sup>0</sup> and *ω* denote the amplitude and frequency of the dynamic load, respectively. The effective material properties for CNTRC beams are given below [27,57]:

$$E\_{11} = \eta\_1 V\_{\text{CNT}} E\_{11}^{\text{CNT}} + V\_p E^p \tag{2}$$

$$\frac{\eta\_{2}}{E\_{22}} = \frac{V\_{CNT}}{E\_{22}^{CNT}} + \frac{V\_{p}}{E^{p}} \tag{3}$$

$$\frac{\eta\_3}{G\_{12}} = \frac{V\_{\rm CNT}}{G\_{12}^{\rm CNT}} + \frac{V\_p}{G^p} \tag{4}$$

$$V\_{\rm CNT} + V\_p = 1\tag{5}$$

$$v = V\_{\text{CNT}} v^{\text{CNT}} + V\_p v^p \tag{6}$$

$$
\rho = V\_{\text{CNT}} \rho^{\text{CNT}} + V\_p \rho^p \tag{7}
$$

where *E*, *G*, *v*, and *ρ* are the material properties that represent the Young's modulus of elasticity, shear modulus, Poisson's ratio, and density, respectively. The superscripts of CNT and p respectively symbolize the related material properties of the carbon nanotube and polymer matrix. *η*1, *η*2, *η*<sup>3</sup> can be indicated the efficiency parameters of CNT. Additionally, *VCNT* and *Vp* define the volume fractions for CNT and polymer matrix, respectively.

The axial strain (*εz*) and shear strain (*γzy*) are given according to the first-order shear deformation beam theory, as follows:

$$
\varepsilon\_z = \frac{\partial u\_0}{\partial z} - \Upsilon \frac{\partial \mathcal{Q}}{\partial z} \tag{8a}
$$

$$
\gamma\_{zy} = \frac{\partial u\_0}{\partial y} + \frac{\partial v\_0}{\partial z} \tag{8b}
$$

where *<sup>u</sup>*0, *<sup>v</sup>*0, and ∅ are axial and vertical displacements, and rotation, respectively. The constitute relation is given below:

$$
\sigma\_z = E(Y) \left[ \frac{\partial u\_0}{\partial z} - Y \frac{\partial \mathcal{D}}{\partial z} \right] \tag{9a}
$$

$$
\sigma\_{zy} = G(Y) K\_s \left[ \frac{\partial v\_0}{\partial z} - \mathcal{Q} \right] \tag{9b}
$$

where *E* is the Young's modulus, *G* is the shear modulus, *σ<sup>z</sup>* is the normal stress, *σzy* is the shear stress, and *Ks* is the shear correction factor.

The strain energy (*U*i), the kinetic energy (*K*), the dissipation function, and potential energy of the external loads (*U*e) are presented as follows:

$$dI\_{i} = \frac{1}{2} \int\_{0}^{L} \left[ A\_{0} \left( \frac{\partial u\_{0}}{\partial z} \right)^{2} - 2A\_{1} \frac{\partial u\_{0}}{\partial z} \frac{\partial \mathcal{Q}}{\partial z} + A\_{2} \left( \frac{\partial \mathcal{Q}}{\partial z} \right)^{2} \right] dZ + \frac{1}{2} \int\_{0}^{L} K\_{i} B\_{0} \left[ \left( \frac{\partial v\_{0}}{\partial z} \right)^{2} - 2 \frac{\partial v\_{0}}{\partial z} \mathcal{Q} + \mathcal{B}^{2} \right] dZ \tag{10a}$$

$$K = \frac{1}{2} \int\_0^L \left( I\_0 \left( \frac{\partial u\_0}{\partial t} \right)^2 - 2I\_1 \left( \frac{\partial u\_0}{\partial t} \right) \left( \frac{\partial \mathcal{Q}}{\partial t} \right) + I\_2 \left( \frac{\partial \mathcal{Q}}{\partial t} \right)^2 + I\_0 \left( \frac{\partial v\_0}{\partial t} \right)^2 \right) dZ \tag{10b}$$

$$\mathcal{U}\_{\varepsilon} = -\mathcal{Q}(t) \, v\left(z\_{p'}t\right) \tag{10c}$$

$$E(A\_{0\prime}, A\_{1\prime}A\_2) = \int\_A E(Y) \{1, Y, Y^2\} dA\_{\prime} \tag{11a}$$

$$B\_0 = \int\_A G(Y)dA\tag{11b}$$

$$\rho(I\_0, I\_1, I\_2) = \int\_A \rho(\mathbf{Y}) \Big( \mathbf{1}, \mathbf{Y}, \mathbf{Y}^2 \Big) dA \tag{11c}$$

The Lagrangian functional of the problem is presented as:

$$I = K - \left(\mathcal{U}\_i + \mathcal{U}\_\ell\right) \tag{12}$$

In the solution of the problem in the Ritz method, the approximate solution is given as a series of *i* terms, as in the following:

$$u\_0(z,t) = \sum\_{i=1}^{\infty} a\_i\left(t\right)a\_i(z) \tag{13a}$$

$$v\_0(z,t) = \sum\_{i=1}^{\infty} b\_i(t)\beta\_i(z) \tag{13b}$$

$$\varpi(z,t) = \sum\_{i=1}^{\infty} c\_i(t)\gamma\_i(z) \tag{13c}$$

where *ai*, *bi*, and *ci* are the unknown coefficients, and *αi*(*z*, *t*), *βi*(*z*, *t*), and *γi*(*z*, *t*) are the coordinate functions depending on the end conditions over the interval [*0*, *L*]. The coordinate functions for the simply supported beam are given as algebraic polynomials:

$$\mathfrak{a}\_{i}(z) = z^{i} \tag{14a}$$

$$\beta\_{\dot{l}}(z) = (L - z)z^{\dot{l}} \tag{14b}$$

$$
\gamma\_i(z) = z^{(i-1)} \tag{14c}
$$

where *i* indicates the number of polynomials involved in the admissible functions. After substituting Equation (7) into Equation (4), and then using the Lagrange's equation, the following equation can be derived:

$$
\frac{\partial \dot{q}\_i}{\partial q\_i} - \frac{\partial t}{\partial t} \frac{\partial \dot{q}\_i}{\partial \dot{q}\_i} = 0 \tag{15}
$$

where *qi* is the unknown coefficients which are *ai*, *bi*, and *ci*. After implementing the Lagrange procedure, the motion equation of the problem was obtained as follows:

$$\left\{ \left[ K \right] \left\{ q(t) \right\} + \left[ M \right] \left\{ \ddot{q}(t) \right\} \right\} = \left\{ F(t) \right\} \tag{16}$$

where [*K*], [*M*], and {*F*(*t*)} are the stiffness matrix, mass matrix, and load vector, respectively. The details of these expressions are given as follows:

$$[K] = \begin{bmatrix} K\_{11} & K\_{12} & K\_{13} \\ K\_{21} & K\_{22} & K\_{23} \\ K\_{31} & K\_{32} & K\_{33} \end{bmatrix} \tag{17}$$

$$K\_{ij}^{11} = \sum\_{i=1}^{n} \sum\_{j=1}^{n} \int\_{0}^{L} A\_0 \frac{\partial \alpha\_i}{\partial z} \frac{\partial \alpha\_j}{\partial z} dz \tag{18a}$$

$$K\_{ij}^{12} = 0,\tag{18b}$$

$$K\_{ij}^{13} = -\sum\_{i=1}^{n} \sum\_{j=1}^{n} \int\_{0}^{L} A\_{1} \frac{\partial \alpha\_{i}}{\partial z} \frac{\partial \gamma\_{j}}{\partial z} dz \tag{18c}$$

$$K\_{ij}^{21} = 0,\tag{18d}$$

$$K\_{ij}^{22} = \sum\_{i=1}^{n} \sum\_{j=1}^{n} \int\_{0}^{L} K\_s B\_0 \frac{\partial \beta\_i}{\partial z} \frac{\partial \beta\_j}{\partial z} dz \tag{18e}$$

$$K\_{ij}^{23} = -\sum\_{i=1}^{n} \sum\_{j=1}^{n} \int\_{0}^{L} K\_{s} B\_{0} \frac{\partial \beta\_{i}}{\partial z} \gamma\_{j} dz \tag{18f}$$

$$K\_{ij}^{31} = -\sum\_{i=1}^{n} \sum\_{j=1}^{n} \int\_{0}^{L} A\_{1} \frac{\partial \gamma\_{i}}{\partial z} \frac{\partial \alpha\_{j}}{\partial z} dz \tag{18g}$$

$$K\_{ij}^{32} = -\sum\_{i=1}^{n} \sum\_{j=1}^{n} \int\_{0}^{L} K\_{s} B\_{0} \gamma\_{i} \frac{\partial \beta\_{j}}{\partial z} dz \tag{18h}$$

$$K\_{ij}^{33} = \sum\_{i=1}^{n} \sum\_{j=1}^{n} \int\_{0}^{L} A\_2 \frac{\partial \gamma\_i}{\partial z} \frac{\partial \gamma\_j}{\partial z} + \sum\_{i=1}^{n} \sum\_{j=1}^{n} \int\_{0}^{L} K\_s B\_0 \gamma\_i \gamma\_j dz \tag{18i}$$

$$\begin{aligned} [M] = \begin{bmatrix} M\_{11} & M\_{12} & M\_{13} \\ M\_{21} & M\_{22} & M\_{23} \\ M\_{31} & M\_{32} & M\_{33} \end{bmatrix} \end{aligned} \tag{19}$$

where

$$M\_{11} = \sum\_{i=1}^{n} \sum\_{j=1}^{n} \int\_{0}^{L} I\_0 \alpha\_i \alpha\_j dz \tag{20a}$$

$$M\_{12} = 0\tag{20b}$$

$$M\_{13} = -\sum\_{i=1}^{n} \sum\_{j=1}^{n} \int\_{0}^{L} I\_1 a\_i \gamma\_j dz \tag{20c}$$

$$M\_{21} = 0\tag{20d}$$

$$M\_{22} = \sum\_{i=1}^{n} \sum\_{j=1}^{n} \int\_{0}^{L} I\_0 \beta\_i \beta\_j dz \tag{20e}$$

$$M\_{2\overline{3}} = 0\tag{20f}$$

$$M\_{31} = -\sum\_{i=1}^{n} \sum\_{j=1}^{n} \int\_{0}^{L} I\_1 \gamma\_i \alpha\_j dz \tag{20g}$$

$$M\_{32} = 0\tag{20h}$$

$$M\_{33} = \sum\_{i=1}^{n} \sum\_{j=1}^{n} \int\_{0}^{L} I\_2 \gamma\_i \gamma\_j dz \tag{20i}$$

$$\{F(t)\} = Q\beta\_{\dot{j}}\tag{21}$$

The constitutive equation of motions was solved by implementing implicit Newmark average acceleration method with *α* = 0.5 and *β* = 0.25 in the time domain. By this procedure, the dynamic problem will be transferred to a system of static problems in each step, as in the following:

$$\left[ \overline{K}(t, X) \right] \{ d\_n \}\_{j+1} = \left\{ \overline{F}(t) \right\} \tag{22}$$

in which

$$\left[\overline{\mathbf{K}}(t, \mathbf{X})\right] = \left[\mathbf{K}\right] + \frac{\left[M\right]}{\beta \Delta t^2} + \frac{\left[\mathbf{C}\right] \alpha}{\beta \Delta t} \tag{23a}$$

$$\left\{\bar{F}(t)\right\} = \left\{F(t)\right\}\_{j+1} + B\_1 \left\{d\_n\right\}\_j + B\_2 \left\{\dot{d}\_n\right\}\_j + B\_3 \left\{\ddot{d}\_n\right\}\_j \tag{23b}$$

$$B\_1 = \frac{[M]}{\beta \Delta t^2}, \ B\_2 = \frac{[M]}{\beta \Delta t}, \ B\_3 = [M] \left(\frac{1}{2\beta} - 1\right) \tag{24}$$

After evaluating {*dn*}*j*+<sup>1</sup> at a time *tj*+<sup>1</sup> = *tj* + Δ*t*, the acceleration and velocity vectors can be expressed as:

$$\left\{\bar{d}\_{n}\right\}\_{j+1} = \frac{1}{\beta \Delta t^{2}} \left( \left\{d\_{n}\right\}\_{j+1} - \left\{d\_{n}\right\}\_{j} \right) - \frac{[M]}{\beta \Delta t} \left\{\bar{d}\_{n}\right\}\_{j} - \left(\frac{a}{2\beta} - 1\right) \left\{\bar{d}\_{n}\right\}\_{j} \tag{25a}$$

$$\left\{\dot{d}\_{n}\right\}\_{j+1} = \left\{\dot{d}\_{n}\right\}\_{j} + \Delta t \left(1 - a\right) \left\{\ddot{d}\_{n}\right\}\_{j} + \Delta t \left\{\ddot{d}\_{n}\right\}\_{j+1} \tag{25b}$$

#### **3. Numerical Results**

In the numerical study, the effects of carbon nanotube volume fraction, aspect ratio, and dynamic parameters on the forced vibration response of CNTRC beams are presented and discussed. The five-point Gauss rule was employed to calculate the integration. In the numerical results, the number of terms is taken as 10. Volume fractions of CNTs as functions of thickness direction for different patterns of CNTs [27] are presented in Table 1. In this table, *V*∗ *CNT* is the given volume fraction of CNTs.

**Table 1.** Volume fractions of CNTs as a function of thickness direction for different distributions of CNTs [27].


In this study, it is notable that CNTs are parallel to the longitudinal direction of the composite beam. Additionally, the efficiency parameters of CNTs for three different values of *V*∗ *CNT* were considered as [23]:

$$
\eta\_1 = 1.2833, \; \eta\_2 = \eta\_3 = 1.0556 \text{ for } V\_{CNT}^\* = 0.12 \tag{26a}
$$

$$
\eta\_1 = 1.3414, \ \eta\_2 = \eta\_3 = 1.7101 \text{ for } V\_{CNT}^\* = 0.17 \tag{26b}
$$

$$
\eta\_1 = 1.3238, \; \eta\_2 = \eta\_3 = 1.7380 \text{ for } V\_{CNT}^\* = 0.28 \tag{26c}
$$

In the numerical results, the following dimensionless displacement was used:

$$
\overline{v} = \frac{E^p b h^3}{PL^3} v
\tag{27}
$$

In the present analysis, the material properties for reinforcement and matrix constituents were [23,27]: *ECNT* <sup>11</sup> = <sup>600</sup> *GPa*, *<sup>E</sup>CNT* <sup>22</sup> = <sup>10</sup> *GPa*, *<sup>G</sup>CNT* <sup>12</sup> = 17.2 *GPa*, *<sup>v</sup>CNT* = 0.19, *ρCNT* = 1400 *kg*/*m*3, *E<sup>p</sup>* = 2.5 *GPa*, *v<sup>p</sup>* = 0.30, and *ρ<sup>p</sup>* = 1190 *kg*/*m*3.

In order to validate the present formulations and analyses, some comparative results are listed in Tables 2 and 3. Firstly, a comparison of non-dimensional fundamental frequencies (λ = *<sup>L</sup>*<sup>2</sup> *<sup>h</sup>* / /*ρa*/*Ea*) of simply supported functionally graded CNT/Aluminum (Al)-alloy composite beams with ANSYS results [19] is presented in Table 2. Here, *k* is the power-law index, and *Ea* and *ρ<sup>a</sup>* represent the elastic modulus and density of pure Al-alloy material, respectively. It can be observed (according to Table 2) that the present results agree well with the ANSYS results [19]. The dimensionless fundamental frequencies (*ω<sup>b</sup>* = *ωL*/ <sup>√</sup>*I*0/*A*0) of simply supported CNTRC beams were calculated with different volume fractions of CNTs for *L/h* = 15 and *VCNT* = 0.12 compared with those of Wattanasakulpong and Ungbhakorn [27], corresponding to the first-order shear deformation theory. To obtain the vibration frequency in this study, the eigenvalue process is implemented in Equation (16). It is seen from Table 3 that the present results are in good agreement with that of the results of Wattanasakulpong and Ungbhakorn [27].

**Table 2.** Comparison of non-dimensional fundamental frequencies of simply supported functionally graded CNT/Al-alloy composite beams with ANSYS results.


**Table 3.** Comparative results for dimensionless fundamental frequencies of a simply supported CNTRC beam for *L/h* = 15, *VCNT* = 0.12.


In order to investigate the effects and compare different reinforcement patterns on dynamic responses, time responses of the simply supported CNTRC beams are presented in Figures 2–4 for volume fractions of CNTs of *VCNT* = 0.12, *VCNT* = 0.17, and *VCNT* = 0.28, respectively. In these figures, the dimensionless midpoint displacements (*vm*) of the beam are obtained within time history for aspect ratio *L/h* = 7 and the external load frequency *ω* = 10 *rd*/*s*. In addition, the dynamical dimensionless displacements of the midpoint (*vm*) and the frequency of the dynamic load (*w*) relations are presented for different reinforcement patterns for *L/h* = 10 and *t* = 1 s in Figures 5–7 for volume fractions of CNTs *VCNT* = 0.12, *VCNT* = 0.17, and *VCNT* = 0.28.

**Figure 2.** Time responses of the CNTRC beam with different reinforcement patterns for *VCNT* = 0.12.

**Figure 3.** Time responses of the CNTRC beam with different reinforcement patterns for *VCNT* = 0.17.

**Figure 4.** Time −responses of the CNTRC beam with different reinforcement pattern for *VCNT* = 0.28.

It is shown, according to Figures 2–10, that the dynamical displacements of the Λ-Beam are the biggest of all. In Λ-Beam, the reinforcements spread only at the bottom surface, not upper surface. However, the reinforcements spread at both surfaces on the UD-Beam and X-Beam. It was known that the upper and lower surfaces of the beam have high stresses and strains. Therefore, the Λ-Beam model had the lowest rigidity in all models. As a result, more displacements occured in the Λ-Beam model. This situation could be observed in Table 3. The vibration frequency of the Λ-Beam was lower than the frequency of the other models. Additionally, the dynamical displacements of the UD-Beam were bigger than those of the X-beam. This is because of the reinforcements spread at both surfaces in the UD- and X-Beams, the X-beam has the biggest specific strength in all patterns. Therefore, dynamic response of the X-Beam is lower than all.

Influence of volume fractions of CNTs on the resonance frequencies of the reinforced composite beam with different distribution patterns are revealed in Table 4. It was found that an increase in volume fractions of CNTs gives rise to an increment in resonance frequencies. Additionally, the highest resonance frequencies occur in the X-Beam, while the lowest ones occur in Λ-Beam.

**Table 4.** The resonance frequencies of reinforced composite beams for various volume fractions of CNTs (*L* = 2 m, *b* = *h* = 0.1 m).


In Figures 5–7, the resonance phenomenon can be observed in the vertical asymptote regions. In the Λ-Beam, the resonance frequency is the lowest for all reinforcement distribution models, because the rigidity of the Λ-Beam is lowest for all. Increasing the volume fractions of CNTs (*VCNT*) yields increased resonance frequency and decreased displacements. It can be interpreted that, by increasing the volume fractions of CNTs, the beam gets more strength. Therefore, the resonance frequencies increase and dynamically displacements decrease naturally.

**Figure 5.** The relationship between the displacements and the frequency of the dynamic load (*w*) for *VCNT* = 0.12.

**Figure 6.** The relationship between the displacements and the frequency of the dynamic load (*w*) for *VCNT* = 0.17.

**Figure 7.** The relationship between the displacements and the frequency of the dynamic load (*w*) for *VCNT* = 0.28.

Figures 8–10 display the frequency of the dynamic load (*w*)-dimensionless vertical displacements relationship for *L/h* = 10 and *t* = 1 *s* for different values of *VCNT* for the UD-beam, Λ-Beam, and X-Beam. It was observed that the increase in values of *VCNT* cause a decrease in the amplitudes of displacements. In the X-beam, the resonance frequencies obtained were the lowest values, in contrast with other values of reinforcement patterns. Another result is that the difference among values of *VCNT* is the highest in the Λ-Beam. It can be concluded that the effects of volume fractions of CNTs were more effective in

Λ-Beams. It shows that the distribution of the reinforcement plays an important role on dynamic responses of CNTRC beams.

**Figure 8.** The relationship between the displacements and the frequency of the dynamic load (*w*) in the simply supported uniformly distributed (UD)-beam for different values of *VCNT*.

**Figure 9.** The relationship between the displacements and the frequency of the dynamic load (*w*) in the simply supported Λ-Beam for different values of *VCNT*.

**Figure 10.** The relationship between of the displacements and the frequency of the dynamic load (*w*) in the simply supported X-Beam for different values of *VCNT*.

#### **4. Conclusions**

The forced vibration response of a simply supported CNTRC beam subjected to a harmonic point load was investigated. It was considered that the composite beam was composed of a polymeric matrix (Poly methyl methacrylate) and a reinforcing material (single-walled carbon nanotubes). Timoshenko beam theory was employed in order to take into consideration the effects of shear deformation. The Ritz and Newmark average acceleration methods were used to obtain the numerical results. The effects of volume fraction and distribution patterns of CNTs, aspect ratio, and dynamic parameters on the forced vibration behavior of CNTRC beam were investigated in detail. It was observed that the greatest dynamical displacements occurred in the Λ-Beam, dependent on having the smallest rigidity. Additionally, it was found that the lowest resonance frequencies were obtained in the X-Beam. In addition, it was revealed that an increase in the values of VCNT gave rise to a decrement in the amplitudes of displacements. Moreover, it was emphasized that the distribution pattern of the reinforcement plays an important role on dynamic responses of CNTRC beams.

**Author Contributions:** Conceptualization, Ö.C. and ¸S.D.A.; methodology, ¸S.D.A.; software, ¸S.D.A.; validation, ¸S.D.A. and B.A.; formal analysis, ¸S.D.A.; investigation, ¸S.D.A.; resources, B.A. and S.D.; data curation, S.D.; writing—original draft preparation, ¸S.D.A. and B.A.; writing—review and editing, Ö.C. and S.D.; visualization, B.A.; supervision, Ö.C.; project administration, Ö.C. and S.D.; funding acquisition, Ö.C. and S.D. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research received no external funding.

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** Data sharing not applicable.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


MDPI St. Alban-Anlage 66 4052 Basel Switzerland Tel. +41 61 683 77 34 Fax +41 61 302 89 18 www.mdpi.com

*Nanomaterials* Editorial Office E-mail: nanomaterials@mdpi.com www.mdpi.com/journal/nanomaterials

MDPI St. Alban-Anlage 66 4052 Basel Switzerland

Tel: +41 61 683 77 34

www.mdpi.com

ISBN 978-3-0365-5777-9