Nonlinear dynamic analysis of functionally graded carbon nanotube-reinforced composite plates using MISQ20 element

Abstract

The main objective of this study is to further extend the mixed integration smoothed quadrilateral element with 20 unknowns of displacement (MISQ20) to investigate the nonlinear dynamic responses of functionally graded carbon nanotube-reinforced composite (FG-CNTRC) plates with four types of carbon nanotube distributions. The smooth finite element method is used to enhance the accuracy of the Q4 element and avoid shear locking without using any shear correction factors. This method yields accurate results even if the element exhibits a concave quadrilateral shape and reduces the error when the element meshing is rough. Additionally, the element stiffness matrix is established by integrating the boundary of the smoothing domains. The motion equation of the FG-CNTRC plates is solved by adapting the Newmark method combined with the Newton–Raphson algorithm. Subsequently, the calculation program is coded in the MATLAB software and verified by comparing it with other published solutions. Finally, the effects of the input parameters on the nonlinear vibration of the plates are investigated.

Graphical abstract

Keywords

carbon nanotube / MISQ20 / FG-CNTRC plate / nonlinear vibration / nonlinear dynamic analysis / SFEM

1 Introduction

Currently, carbon nanotubes (CNTs) are widely employed in engineering and have piqued the interest of scientists worldwide. Thostenson et al. [1] improved the characteristics of composite materials using specific methods. CNTs are considered promising reinforcements for polymer composites [2]. Wagner et al. [3,4] experimentally measured the interfacial strength CNT–polymer by performing numerous pull-out tests involving single CNTs and discovered the effectiveness of CNTs in improving the polymer. Gou et al. [5] simulated the effect of chemical links on the shear modulus of CNTs, whereas Frankland et al. [6] numerically and experimentally investigated the interfacial bonding of single-walled carbon nanotube-reinforced composite (CNTRC). Ma et al. [7] investigated a method to enhance amino-functionalized CNTs in a material matrix.

The mechanical properties of CNTs have been investigated using several methods, such as analytical methods, experimental solutions, and numerical approaches. Coleman et al. [8] provided a general review of the mechanical behavior of functionally graded carbon nanotube-reinforced composite (FG-CNTRC) plates. The authors conducted experiments and discovered that adding small amounts of CNT to the original material significantly increased the stiffness of the plates [9,10]. Odegard et al. [11,12] performed an equivalent continuum modeling and constitutive modeling of CNTRC polymers via molecular dynamic (MD) simulations. Liu and Chen [13] and Hu et al. [14] used a three-dimensional finite element method (FEM) to consider the macroscopic elastic properties of CNTs. Frankland et al. [15] and Griebel and Hamaekers [16] performed MD simulations to calculate the different elastic moduli of CNTs. Wernik and Meguid [17] investigated the nonlinear behavior of CNT–polymers. CNTs are advanced materials used in engineering applications, including CNTRC structures; therefore, the mechanical behavior of the latter must be investigated. Wuite and Adali [18] used a length-scale method to compute the static response of a nanocomposite beam. Vodenitcharova and Zhang [19] modified a continuum model to analyze the mechanical behavior of nanocomposite beams reinforced with single-walled carbon nanotubes (SWCNTs). Ray and Batra [20] controlled the vibration of SWCNT structures using piezoelectric materials. Formica et al. [21] applied the Tanaka theory to compute the vibrations of CNTRC plates. They confirmed that the vibrations can be customized via a specified external excitation. Arani et al. [22] examined the buckling of plates using both an exact solution and a numerical approach. Shen [23] investigated the nonlinear statics of plates using HSDT. Wang and Shen [24] computed the large-magnitude oscillations of plates lying on an elastic substrate. Furthermore, Wang and Shen [25] computed the nonlinear mechanical behavior of sandwich plates using an analytical method.

Using the FEM to analyze FG-CNTRC structures allows us to summarize a few previous studies, such as those of Zhang et al. [26] employed the element-free improved moving least-squares Ritz method to examine the vibration of FG-CNTRC plates. Natarajan et al. [27] employed the QUAD-8 element combined with HSDT to examine the static and vibrational properties of FG-CNTRC plates. Using the QUAD-8 element, Sankar et al. analyzed the flutter of sandwich plates [28] and doubly curved sandwich panels [29] as well as computed the nonlinear buckling of spherical/conical shells [30]. Rodrigues et al. [31] examined the static bending of plates using the MITC4 element.

To improve computational efficiency, Hughes et al. [32] and Borden et al. [33] introduced isogeometric analysis (IGA), which is based on computer-aided design (CAD). In recent years, the use of IGA for structural analysis has garnered the attention of many scientists. Phung-Van et al. [34] performed IGA to investigate the bending and free vibration of porous functionally graded (FG) nanoplates. Thanh et al. [35] combined IGA with modified couple stress theory to examine the static of FG microplates. Furthermore, they performed IGA combined with a nonlocal strain gradient to analyze the mechanical behavior of FG sandwiched nanoplates [36]. Moreover, some results of structural analysis obtained via IGA have been reported [37,38]. To overcome the limitations of FEMs with arbitrary structural shapes, Guo and Zheng improved the numerical manifold method to analyze shells [39] and Kirchhoff’s plates with arbitrary geometries [40]. Additionally, they used deep neural networks [41], applied machine learning [42], and provided a MATLAB code [43] to avoid classical discretization. Their approach can be adopted in future structural analyses.

Ganapathi et al. [44] investigated the nonlinear vibration of structures. Kant and Kommineni [45] used a refined plate formulation to examine the vibrations of composite structures. Dash [46] calculated the nonlinear vibration of sandwich structures using a first-order FEM. Ather [47] and Parhi et al. [48] calculated the nonlinear free vibration of composite structures using an FEM. Additionally, flat elements have been widely used [49–54], which can be easily combined with other elements to improve computational performance. Results show that flat elements are convenient for solving geometrically nonlinear issues where the structural behavior based on the iterative solution must be addressed and saved with numerous history variables. Other solutions for the nonlinear behavior analysis of structures are available in the literature [55–57].

Scientists worldwide have investigated the mechanical behavior of FG-CNTRC structures using different methods, including analytical, numerical, and meshless methods, as well as IGA. Whereas linear problems were typically prioritized, nonlinear problems were rarely mentioned and analytical methods were often used. In addition, using the Q4 element facilitates numerical computation and yields the initial results for analyzing the mechanical behavior of structures. Therefore, we were motivated to use the mixed integration smoothed quadrilateral element with 20 unknowns of displacement (MISQ20) element based on the Q4 element to establish a small strain–large deformation formulation. Subsequently, we computed the nonlinear vibration of the plates using the Newmark method combined with the Newton–Raphson iteration. Based on the smooth finite element method, the stiffness matrices were obtained by integrating the boundary of the smoothing domains. This integration technique provides high accuracy even if distorted elements or coarse meshes are used and reduces the operating time of the calculation program. Some examples are provided herein to show the nonlinear vibration results of the plates.

2 CNT-reinforced composites

FG-CNTRC plates with four types of CNTs were investigated in this study, as shown in Fig.1. Using the Mori–Tanaka method, the mechanical properties of FG-CNTRC plates were obtained by combining straight CNTs and an isotropic polymer as follows [23]:

(1)\{\begin{array}{rl}& E_{11}={\eta }_1{V}_{\mathrm{C}\mathrm{N}\mathrm{T}}{E}_{11}^{\mathrm{C}\mathrm{N}\mathrm{T}}+V_{\mathrm{m}}{E}^{\mathrm{m}},\\ & \frac{{\eta }_2}{E_{22}}=\frac{V_{\mathrm{C}\mathrm{N}\mathrm{T}}}{E_{22}^{\mathrm{C}\mathrm{N}\mathrm{T}}}+\frac{V_{\mathrm{m}}}{E^{\mathrm{m}}},\\ & \frac{{\eta }_3}{G_{12}}=\frac{V_{\mathrm{C}\mathrm{N}\mathrm{T}}}{G_{22}^{\mathrm{C}\mathrm{N}\mathrm{T}}}+\frac{V_{\mathrm{m}}}{G^{\mathrm{m}}},\\ & \rho ={\rho }_{\mathrm{C}\mathrm{N}\mathrm{T}}{V}_{\mathrm{C}\mathrm{N}\mathrm{T}}+{\rho }_{\mathrm{m}}{V}_{\mathrm{m}},\end{array}

where symbols “CNT” and “m” represent the properties of the CNTs and original polymer, respectively. The elastic modulus, shear modulus, and material density of the isotropic matrix are denoted as E^{\mathrm{m}}, G^{\mathrm{m}}, and \rho, respectively. In addition, the efficiency parameters {\eta }_1, {\eta }_2, and {\eta }_3 are provided in Ref. [57].

In particular, {\eta }_1=0.137,{\eta }_2=1.022,{\eta }_3=0.715 for V_{\mathrm{C}\mathrm{N}\mathrm{T}}^{\ast }=0.12 (12%); {\eta }_1=0.142,{\eta }_2=1.626,{\eta }_3=1.138 for V_{\mathrm{C}\mathrm{N}\mathrm{T}}^{\ast }=0.17 (17%), and {\eta }_1=0.141,{\eta }_2=1.585,{\eta }_3= 1.109 for V_{\mathrm{C}\mathrm{N}\mathrm{T}}^{\ast }=0.28 (28%). V_{\mathrm{C}\mathrm{N}\mathrm{T}} and V_{\mathrm{m}} denote the volume fractions of the CNT and polymer, respectively, and the relationship between V_{\mathrm{C}\mathrm{N}\mathrm{T}} and V_{\mathrm{m}} is expressed as follows [58]:

(2)V_{\mathrm{C}\mathrm{N}\mathrm{T}}+V_{\mathrm{m}}=1.

The Poisson’s ratio is expressed by

(3)v_{12}=V_{\mathrm{C}\mathrm{N}\mathrm{T}}^{\ast }{v}_{12}^{\mathrm{C}\mathrm{N}\mathrm{T}}+V_{\mathrm{m}}^{}{v}_{}^{\mathrm{m}}.

The volume of the CNT is estimated as follows:

(4)V_{\mathrm{C}\mathrm{N}\mathrm{T}}=\{\begin{array}{rlr}& V_{\mathrm{C}\mathrm{N}\mathrm{T}}^{\ast },& \text{UD,}\\ & \left(1+\frac{2\mathit{\text{z}}}{h}\right)V_{\mathrm{C}\mathrm{N}\mathrm{T}}^{\ast },& \text{FG-V,}\\ & 2\left(1-\frac{2\left|\mathit{\text{z}}\right|}{h}\right)V_{\mathrm{C}\mathrm{N}\mathrm{T}}^{\ast },& \text{FG-O,}\\ & \frac{4\left|\mathit{\text{z}}\right|}{h}{V}_{\mathrm{C}\mathrm{N}\mathrm{T}}^{\ast },& \text{FG-X,}\end{array}

with

(5)V_{\mathrm{C}\mathrm{N}\mathrm{T}}^{\ast }=\frac{{\beta }_{\mathrm{C}\mathrm{N}\mathrm{T}}}{{\beta }_{\mathrm{C}\mathrm{N}\mathrm{T}}+\left({\rho }_{\mathrm{C}\mathrm{N}\mathrm{T}}/{\rho }_{\mathrm{m}}\right)\left(1-{\beta }_{\mathrm{C}\mathrm{N}\mathrm{T}}\right)},

where {\beta }_{\mathrm{C}\mathrm{N}\mathrm{T}} denotes the mass fraction of the CNT; and {\rho }_{\mathrm{m}} and {\rho }_{\mathrm{C}\mathrm{N}\mathrm{T}} represent the densities of the polymer and CNT, respectively.

3 Finite element formulations

Based on FSDT, the displacement field of plates is expressed as [59]

(6)\mathit{q}=\{\begin{array}{c}{u}_1(x,y,\mathit{\text{z}})=u_0(x,y,0)+\mathit{\text{z}}{\theta }_x(x,y,0),\\ u_2(x,y,\mathit{\text{z}})=v_0(x,y,0)+\mathit{\text{z}}{\theta }_y(x,y,0),\\ u_3(x,y,\mathit{\text{z}})=w_0(x,y,0),\end{array}

where u_0,v_0,w_0,{\theta }_x, and {\theta }_y are unknown displacements.

The stress–strain relationship is expressed as

(7)\{\begin{array}{c}{\epsilon }_x=u_{1,x}=u_{0,x}+\mathit{\text{z}}{\theta }_{x,x}+\frac{1}{2}{w}_{0,x}^2,\\ {\epsilon }_y=u_{2,y}=v_{0,y}+\mathit{\text{z}}{\theta }_{y,y}+\frac{1}{2}{w}_{0,y}^2,\\ {\gamma }_{xy}=u_{1,y}+u_{2,x}=u_{0,y}+v_{0,x}+w_{0,x}{w}_{0,y},\\ {\gamma }_{x\mathit{\text{z}}}=w_{0,x}+u_{1,\mathit{\text{z}}}=w_{0,x}+{\theta }_x,\\ {\gamma }_{y\mathit{\text{z}}}=w_{0,y}+u_{2,\mathit{\text{z}}}=w_{0,y}+{\theta }_y,\end{array}

By substituting Eq. (6) into Eq. (7), the in-plane strain vector can be rewritten as

(8)\mathit{\epsilon }=\left\{\begin{array}{c}{\epsilon }_x\\ {\epsilon }_y\\ {\gamma }_{xy}\end{array}\right\}={\mathit{\epsilon }}_{\mathrm{L}}^0+{\mathit{\epsilon }}_{\mathrm{N}}^0+\mathit{\text{z}}\mathit{\kappa },\mathit{\gamma }=\left\{\begin{array}{c}{\gamma }_{y\mathit{\text{z}}}\\ {\gamma }_{\mathit{\text{z}}x}\end{array}\right\},

where

(9){\mathit{\epsilon }}_{\mathrm{L}}^0={\left\{\begin{array}{ccc}{u}_{0,x}& v_{0,y}& u_{0,y}+v_{0,x}\end{array}\right\}}^{\mathrm{T}},{\mathit{\epsilon }}_{\mathrm{N}}^0=\frac{1}{2}{\left\{\begin{array}{ccc}{w}_{0,x}^2& w_{0,y}^2& 2w_{0,x}{w}_{0,y}\end{array}\right\}}^{\mathrm{T}},

(10)\mathit{\kappa }={\left\{\begin{array}{ccc}{\theta }_{x,x}& {\theta }_{y,y}& {\theta }_{x,y}+{\theta }_{y,x}\end{array}\right\}}^{\mathrm{T}},\mathit{\gamma }={\left\{\begin{array}{cc}{w}_{0,x}+{\theta }_x& w_{0,y}+{\theta }_y\end{array}\right\}}^{\mathrm{T}},

Subsequently, the nonlinear membrane strain–displacement vector is rewritten as follows:

(11){\mathit{\epsilon }}_{\mathrm{N}}^0=\frac{1}{2}\underset{\mathit{H}}{\underbrace{\left[\begin{array}{cc}{w}_{0,x}& 0\\ 0& w_{0,y}\\ w_{0,y}& w_{0,x}\end{array}\right]}}\underset{\mathit{Y}}{\underbrace{\left[\begin{array}{c}{w}_{0,x}\\ w_{0,y}\end{array}\right]}}=\frac{1}{2}\mathit{H}\mathit{Y},

where Y is the slope vector.

Based on Hooke’s law, the stress–strain relationship is written as follows:

(12)\left\{\begin{array}{c}{\sigma }_x\\ {\sigma }_y\\ {\tau }_{xy}\end{array}\right\}=\underset{{\mathit{Q}}_{\mathrm{b}}}{\underbrace{\left[\begin{array}{ccc}{q}_{11}& q_{12}& 0\\ q_{21}& q_{22}& 0\\ 0& 0& q_{66}\end{array}\right]}}\left\{\begin{array}{c}{\epsilon }_x\\ {\epsilon }_y\\ {\gamma }_{xy}\end{array}\right\},

where

(13)\left\{\begin{array}{c}{\tau }_{x\mathit{\text{z}}}\\ {\tau }_{y\mathit{\text{z}}}\end{array}\right\}=\underset{{\mathit{Q}}_{\mathrm{s}}}{\underbrace{\left[\begin{array}{cc}{q}_{55}& 0\\ 0& q_{44}\end{array}\right]}}\left\{\begin{array}{c}{\gamma }_{x\mathit{\text{z}}}\\ {\gamma }_{y\mathit{\text{z}}}\end{array}\right\},

where

q_{11}(\mathit{\text{z}})=q_{22}(\mathit{\text{z}})=\frac{E_{11}}{1-v_{12}},q_{12}(\mathit{\text{z}})=v_{21}{q}_{11}(\mathit{\text{z}}),

(14)q_{44}(\mathit{\text{z}})=q_{55}(\mathit{\text{z}})=q_{66}(\mathit{\text{z}})=\frac{E_{11}}{2\left(1+v_{12}\right)}.

A weak form of the plates for dynamic analysis is presented as follows:

(15)\begin{array}{rl}& {\int }_{\mathrm{\Omega }}\delta {\dot{\mathit{q}}}^{\mathrm{T}}\mathit{m}\dot{\mathit{q}}\mathrm{d}\mathrm{\Omega }+{\int }_{\mathrm{\Omega }}\delta {\mathit{\epsilon }}^{\mathrm{T}}{\mathit{D}}_{\mathrm{b}}\mathit{\epsilon }\mathrm{d}\mathrm{\Omega }+{\int }_{\mathrm{\Omega }}\delta {\mathit{\gamma }}^{\mathrm{T}}{\mathit{D}}_{\mathrm{s}}\mathit{\gamma }\mathrm{d}\mathrm{\Omega }\\ -& {\int }_{\mathrm{\Omega }}\delta {\mathit{q}}^{\mathrm{T}}{\mathit{N}}_{\mathrm{w}}p(t)\mathrm{d}\mathrm{\Omega }=0,\end{array}

where

(16)\mathit{\epsilon }=\left[\begin{array}{c}{\mathit{\epsilon }}_0\\ \mathit{\kappa }\end{array}\right]=\left[\begin{array}{c}{\mathit{\epsilon }}_{\mathrm{L}}^0+{\mathit{\epsilon }}_{\mathrm{N}}^0\\ \mathit{\kappa }\end{array}\right],{\mathit{D}}_{\mathrm{b}}=\left[\begin{array}{cc}\mathit{A}& \mathit{B}\\ \mathit{B}& \mathit{D}\end{array}\right],{\mathit{D}}_{\mathrm{s}}={\int }_{-h/2}^{h/2}{\mathit{Q}}_{\mathrm{s}}\mathrm{d}\mathit{\text{z}},

with

(17)\left(\mathit{A},\mathit{B},\mathit{D}\right)={\int }_{-h/2}^{h/2}\left(1,\mathit{\text{z}},{\mathit{\text{z}}}^2\right){\mathit{Q}}_{\mathrm{b}}\mathrm{d}\mathit{\text{z}}.

The inertial component matrix is expressed as

(18)\begin{array}{rl}& \mathit{m}=\left[\begin{array}{ccccc}{I}_0& 0& 0& I_1& 0\\ & I_0& 0& 0& I_1\\ & & I_0& 0& 0\\ & & & I_2& 0\\ sym& & & & I_2\end{array}\right],\\ & (I_0,I_1,I_2)={\int }_{-h/2}^{h/2}\rho (1,\mathit{\text{z}},{\mathit{\text{z}}}^2)\mathrm{d}\mathit{\text{z}}.\end{array}

Using a four-node plate element with five DOFs per node, the displacement of the elements can be defined as follows:

(19){\mathit{q}}_{\mathrm{e}}={\left\{\begin{array}{cccc}{\mathit{q}}_{\mathrm{e}}^{\text{1}}& {\mathit{q}}_{\mathrm{e}}^2& {\mathit{q}}_{\mathrm{e}}^3& {\mathit{q}}_{\mathrm{e}}^4\end{array}\right\}}^{\mathrm{T}}.

In the above,

(20){\mathit{q}}_{\mathrm{e}}^i={\left\{\begin{array}{ccccc}{u}_0& v_0& w_0& {\theta }_x& {\theta }_y\end{array}\right\}}^{\mathrm{T}},i=1,2,3,4,

where

(21)\{\begin{array}{rl}& u_0=\sum _{i=1}^4{L}_i\left(\xi ,\eta \right)u_{0i},\\ & v_0=\sum _{i=1}^4{L}_i\left(\xi ,\eta \right)v_{0i},\\ & w_0=\sum _{i=1}^4{L}_i\left(\zeta ,\eta \right)w_{0i},\\ & {\theta }_x=\sum _{i=1}^4{L}_i\left(\xi ,\eta \right){\theta }_{xi},\\ & {\theta }_y=\sum _{i=1}^4{L}_i\left(\xi ,\eta \right){\theta }_{yi},\end{array}

where L_i (i = 1,2,3,4) is the shape function [59].

By substituting the displacement component in Eq. (21) in Eq. (9), we obtain

(22)\{\begin{array}{c}{\mathit{\epsilon }}_{\mathrm{L}}^0=\left({\mathit{B}}_{}^0+\mathit{\text{z}}{\mathit{B}}_{}^1\right){\mathit{q}}_{\mathrm{e}},{\mathit{\epsilon }}_{\mathrm{N}}^0={\mathit{B}}^{\mathrm{N}}{\mathit{q}}_{\mathrm{e}},\\ \mathit{\gamma }={\mathit{B}}_{}^2{\mathit{q}}_{\mathrm{e}},\end{array}

where {\mathit{B}}_{}^0,{\mathit{B}}_{}^1, and {\mathit{B}}_{}^2 are expressed as follows:

(23)\begin{array}{c}{\mathit{B}}_{}^0=\left[\begin{array}{cccc}{\mathit{B}}_1^0& {\mathit{B}}_2^0& {\mathit{B}}_3^0& {\mathit{B}}_4^0\end{array}\right],\\ {\mathit{B}}_{}^1=\left[\begin{array}{cccc}{\mathit{B}}_1^1& {\mathit{B}}_2^1& {\mathit{B}}_3^1& {\mathit{B}}_4^1\end{array}\right],\\ {\mathit{B}}_{}^2=\left[\begin{array}{cccc}{\mathit{B}}_1^2& {\mathit{B}}_2^2& {\mathit{B}}_3^2& {\mathit{B}}_4^2\end{array}\right],\end{array}

where

(24)\begin{array}{rl}& {\mathit{B}}_i^0=\left[\begin{array}{ccccc}{L}_{i,x}& 0& 0& 0& 0\\ 0& L_{i,y}& 0& 0& 0\\ L_{i,y}& L_{i,x}& 0& 0& 0\end{array}\right],{\mathit{B}}_i^1=\left[\begin{array}{ccccc}0& 0& 0& L_{i,x}& 0\\ 0& 0& 0& 0& L_{i,y}\\ 0& 0& 0& L_{i,y}& L_{i,x}\end{array}\right],\\ & {\mathit{B}}_i^2=\left[\begin{array}{ccccc}0& 0& L_{i,y}& 0& L_i\\ 0& 0& L_{i,x}& L_i& 0\end{array}\right],\\ & {\mathit{B}}_i^N=\frac{1}{2}\left[\begin{array}{cc}{w}_{0,x}& 0\\ 0& w_{0,y}\\ w_{0,y}& w_{0,x}\end{array}\right]\left[\begin{array}{ccccc}0& 0& L_{i,x}& 0& 0\\ 0& 0& L_{i,y}& 0& 0\end{array}\right]=\frac{1}{2}\mathit{H}\cdot \mathit{G}.\end{array}

By replacing Eqs. (19) and (22) into Eq. (15), the nonlinear vibration equation for the plate element can be written as follows:

(25){\mathit{M}}_{\mathrm{e}}^{}{\ddot{\mathit{q}}}_{\mathrm{e}}+\left({\mathit{K}}_{\mathrm{e}}^{\mathrm{L}}+{\mathit{K}}_{\mathrm{e}}^{\mathrm{N}}\right){\mathit{q}}_{\mathrm{e}}={\mathit{F}}_{\mathrm{e}}.

In the equation above, the element stiffness matrices are

(26)\begin{array}{rl}{\mathit{K}}_{\mathrm{e}}^{\mathrm{L}}=& \underset{\mathrm{\Omega }}{\int }{\left({\mathit{B}}_{}^0\right)}^{\mathrm{T}}{\mathit{A}\mathit{B}}_{}^0\mathrm{d}\mathrm{\Omega }+\underset{\mathrm{\Omega }}{\int }{\left(B_{}^0\right)}^{\mathrm{T}}\mathit{B}{\mathit{B}}^1\mathrm{d}\mathrm{\Omega }+\underset{\mathrm{\Omega }}{\int }{\left({\mathit{B}}_{}^1\right)}^{\mathrm{T}}{\mathit{B}\mathit{B}}_{}^0\mathrm{d}\mathrm{\Omega }\\ & +\underset{\mathrm{\Omega }}{\int }{\left({\mathit{B}}_{}^1\right)}^{\mathrm{T}}{\mathit{D}\mathit{B}}_{}^1\mathrm{d}\mathrm{\Omega }+\underset{\mathrm{\Omega }}{\int }{\left({\mathit{B}}_{}^2\right)}^{\mathrm{T}}{\mathit{D}}_{\mathrm{s}}{\mathit{B}}_{}^2\mathrm{d}\mathrm{\Omega },\end{array}

and

(27)\begin{array}{rl}{\mathit{K}}_{\mathrm{e}}^{\mathrm{N}}=& \underset{\mathrm{\Omega }}{\int }({\left({\mathit{B}}^{\mathrm{N}}\right)}^{\mathrm{T}}\mathit{A}{\mathit{B}}^{\mathrm{N}}+{\left({\mathit{B}}^{\mathrm{N}}\right)}^{\mathrm{T}}\mathit{A}{\mathit{B}}^0+{\left({\mathit{B}}^0\right)}^{\mathrm{T}}\mathit{A}{\mathit{B}}^{\mathrm{N}}\\ & +{\left({\mathit{B}}^{\mathrm{N}}\right)}^{\mathrm{T}}\mathit{B}{\mathit{B}}^1+{\left({\mathit{B}}^{\mathrm{N}}\right)}^{\mathrm{T}}\mathit{B}{\mathit{B}}^1+{\left({\mathit{B}}^{\mathrm{N}}\right)}^{\mathrm{T}}\mathit{A}{\mathit{B}}^1)\mathrm{d}\mathrm{\Omega }.\end{array}

The element mass matrix is written as

(28){\mathit{M}}_{\mathrm{e}}^{}=\underset{\mathrm{\Omega }}{\int }{\mathit{L}}^{\mathrm{T}}\mathit{m}\mathit{L}\mathrm{d}\mathrm{\Omega },

where \mathit{L}=[\begin{array}{cccc}{L}_1{\mathit{I}}_{5\times 5}& L_2{\mathit{I}}_{5\times 5}& L_3{\mathit{I}}_{5\times 5}& L_4{\mathit{I}}_{5\times 5}\end{array}], with {\mathit{I}}_{5\times 5} is the unit matrix of 5 × 5 rank.

The element load vector is expressed as

(29){\mathit{F}}_{\mathrm{e}}^{}=\underset{\mathrm{\Omega }}{\int }{\mathit{L}}^{\mathrm{T}}\mathit{f}\mathrm{d}\mathrm{\Omega },

with \mathit{f}={\left\{\begin{array}{ccccc}0& 0& p(t)& 0& 0\end{array}\right\}}^{\mathrm{T}}.

Considering a quadrilateral element domain, we divided nc into smoothing cells (see Fig.2). The strain field is obtained as follows [58]:

(30){\overline{\mathit{\epsilon }}}_0^{\mathrm{L}}=\frac{1}{A_{\mathrm{C}}}\underset{{\mathrm{\Omega }}_{\mathrm{C}}}{\int }{\mathit{\epsilon }}_0^{\mathrm{L}}(x)\mathrm{d}\mathrm{\Omega },

(31){\overline{\mathit{\epsilon }}}_0^{\mathrm{N}}=\frac{1}{A_{\mathrm{C}}}\underset{{\mathrm{\Omega }}_{\mathrm{C}}}{\int }{\mathit{\epsilon }}_0^{\mathrm{N}}(x)\mathrm{d}\mathrm{\Omega },

(32)\overline{\mathit{\epsilon }}=\frac{1}{A_{\mathrm{C}}}\underset{{\mathrm{\Omega }}_{\mathrm{C}}}{\int }\mathit{\epsilon }(x)\mathrm{d}\mathrm{\Omega },

where A_{\mathrm{C}} is the area of cell {\mathrm{\Omega }}_{\mathrm{C}}.

Based on the divergence principle, the smoothed membrane strain is expressed as [58]

(33)\begin{array}{rl}{\overline{\mathit{\epsilon }}}_0^{\mathrm{L}}& =\frac{1}{A_{\mathrm{C}}}\underset{{\mathrm{\Gamma }}_{\mathrm{C}}}{\int }\mathit{n}(x)\mathit{q}(x)\mathrm{d}{\mathrm{\Gamma }}_{\mathrm{C}}=\frac{1}{A_{\mathrm{C}}}\underset{{\mathrm{\Gamma }}_{\mathrm{C}}}{\int }\sum _{i=1}^4\mathit{n}(x){\mathit{L}}_i(x){\mathit{q}}_{\mathrm{e}}^i\mathrm{d}\mathrm{\Gamma }\\ & =\sum _{i=1}^4\left[\begin{array}{cc}{\overline{\mathit{B}}}_i^0& {\overline{\mathit{B}}}_i^1\end{array}\right]{\mathit{q}}_{\mathrm{e}}^i,\end{array}

where

(34){\overline{\mathit{B}}}_i^0=\frac{1}{A_{\mathrm{C}}}\underset{{\mathrm{\Gamma }}_{\mathrm{C}}}{\int }\left[\begin{array}{ccccc}{L}_i{n}_x& 0& 0& 0& 0\\ 0& L_i{n}_y& 0& 0& 0\\ L_i{n}_y& L_i{n}_x& 0& 0& 0\end{array}\right]\mathrm{d}\mathrm{\Gamma },

(35){\overline{\mathit{B}}}_i^1=\frac{1}{A_{\mathrm{C}}}\underset{{\mathrm{\Gamma }}_{\mathrm{C}}}{\int }\left[\begin{array}{ccccc}0& 0& 0& L_i{n}_x& 0\\ 0& 0& 0& 0& L_i{n}_y\\ 0& 0& 0& L_i{n}_y& L_i{n}_x\end{array}\right]\mathrm{d}\mathrm{\Gamma }.

By performing Gaussian integration, we obtain

(36){\overline{\mathit{B}}}_i^0=\frac{1}{A_{\mathrm{C}}}\sum _{m=1}^4\left[\begin{array}{ccccc}{L}_i{x}_m^{\mathrm{G}}{n}_x& 0& 0& 0& 0\\ 0& L_i{x}_m^{\mathrm{G}}{n}_y& 0& 0& 0\\ L_i{x}_m^{\mathrm{G}}{n}_y& L_i{x}_m^{\mathrm{G}}{n}_x& 0& 0& 0\end{array}\right]l_m^{\mathrm{C}},

(37){\overline{\mathit{B}}}_i^1=\frac{1}{A_{\mathrm{C}}}\sum _{b=1}^4\left[\begin{array}{ccccc}0& 0& 0& L_i(x_b^{\mathrm{G}})n_x& 0\\ 0& 0& 0& 0& L_i(x_b^{\mathrm{G}})n_y\\ 0& 0& 0& L_i(x_b^{\mathrm{G}})n_y& L_i(x_b^{\mathrm{G}})n_x\end{array}\right]l_b^{\mathrm{C}}.

Subsequently, the smoothed nonlinear membrane strain is expressed as

(38){\overline{\mathit{\epsilon }}}_0^{\mathrm{N}}=\sum _{i=1}^4{\overline{\mathit{B}}}_i^{\mathrm{N}}{q}_{\mathrm{e}}^i,

where {\overline{\mathit{B}}}_i^{\mathrm{N}} is calculated as

(39){\mathit{B}}_i^{\mathrm{N}}=\frac{1}{2}\overline{\mathit{H}}{\overline{\mathit{G}}}_i,

in which

(40)\overline{\mathit{H}}={\displaystyle \frac{1}{A_{\mathrm{C}}}}\sum _{i=1}^4\left[\begin{array}{cc}{\displaystyle \sum _{j=1}^4{L}_i{x}_j^{\mathrm{G}}{n}_x{l}_j^{\mathrm{G}}{w}_i}& 0\\ 0& {\displaystyle \sum _{j=1}^4{L}_i{x}_j^{\mathrm{G}}{n}_y{l}_j^{\mathrm{G}}{w}_i}\\ {\displaystyle \sum _{j=1}^4{L}_i{x}_j^{\mathrm{G}}{n}_y{l}_j^{\mathrm{G}}{w}_i}& {\displaystyle \sum _{j=1}^4{L}_i{x}_j^{\mathrm{G}}{n}_x{l}_j^{\mathrm{G}}{w}_i}\end{array}\right],

(41){\overline{\mathit{G}}}_i=\frac{1}{A_{\mathrm{C}}}\sum _{g=1}^4\left[\begin{array}{ccccc}0& 0& L_i{x}_g^{\mathrm{G}}{n}_x& 0& 0\\ 0& 0& L_i{x}_g^{\mathrm{G}}{n}_y& 0& 0\end{array}\right]l_g^{\mathrm{C}}.

The shear strain is approximated in the natural coordinate system as follows:

(42)\left[\begin{array}{c}{\gamma }_{x\mathit{\text{z}}}\\ {\gamma }_{y\mathit{\text{z}}}\end{array}\right]=J^{-1}\left[\begin{array}{cccc}{\displaystyle \frac{1}{2}}(1-\xi )& 0& {\displaystyle \frac{1}{2}}(1+\xi )& 0\\ 0& {\displaystyle \frac{1}{2}}(1-\eta )& 0& {\displaystyle \frac{1}{2}}(1+\eta )\end{array}\right]\left[\begin{array}{c}{\gamma }_{\eta }^{\mathrm{A}}\\ {\gamma }_{\xi }^{\mathrm{B}}\\ {\gamma }_{\eta }^{\mathrm{C}}\\ {\gamma }_{\xi }^{\mathrm{D}}\end{array}\right],

where J is the Jacobian matrix.

The displacement–strain matrix {\overline{\mathit{B}}}_i^2 is approximated as follows:

(43){\overline{\mathit{B}}}_i^2=J^{-1}\left[\begin{array}{ccccc}0& 0& L_{i,\xi }& b_i^{11}{L}_{i,\xi }& b_i^{12}{L}_{i,\xi }\\ 0& 0& L_{i,\eta }& b_i^{21}{L}_{i,\eta }& b_i^{22}{L}_{i,\eta }\end{array}\right].

Finally, the motion equation of the plate element is expressed as

(44){\mathit{M}}_{\mathrm{e}}^{}{\ddot{\mathit{q}}}_{\mathrm{e}}+\left({\overline{\mathit{K}}}_{\mathrm{e}}^{\mathrm{L}}+{\overline{\mathit{K}}}_{\mathrm{e}}^{\mathrm{N}}\right){\mathit{q}}_{\mathrm{e}}={\mathit{F}}_{\mathrm{e}},

with

(45)\begin{array}{rl}{\overline{\mathit{K}}}_{\mathrm{e}}^{\mathrm{L}}=& \underset{\mathrm{\Omega }}{\int }[{\left({\overline{\mathit{B}}}_{}^0\right)}^{\mathrm{T}}{\mathit{A}\overline{\mathit{B}}}_{}^0+{\left({\overline{\mathit{B}}}_{}^0\right)}^{\mathrm{T}}\mathit{B}{\overline{\mathit{B}}}^1+{\left({\overline{\mathit{B}}}_{}^1\right)}^{\mathrm{T}}{\mathit{B}\overline{\mathit{B}}}_{}^0\\ & +{\left({\overline{\mathit{B}}}_{}^1\right)}^{\mathrm{T}}{\mathit{D}\overline{\mathit{B}}}_{}^1+{\left({\overline{\mathit{B}}}_{}^2\right)}^{\mathrm{T}}{\mathit{D}}_{\mathrm{s}}{\overline{\mathit{B}}}_{}^2]\mathrm{d}\mathrm{\Omega },\end{array}

(46)\begin{array}{rl}{\overline{\mathit{K}}}_{\mathrm{e}}^{\mathrm{N}}=& \underset{\mathrm{\Omega }}{\int }[{\left({\overline{\mathit{B}}}^{\mathrm{N}}\right)}^{\mathrm{T}}\mathit{A}{\overline{\mathit{B}}}^N+{\left({\overline{\mathit{B}}}^{\mathrm{N}}\right)}^{\mathrm{T}}\mathit{A}{\overline{\mathit{B}}}^0+{\left({\overline{\mathit{B}}}^0\right)}^{\mathrm{T}}\mathit{A}{\overline{\mathit{B}}}^{\mathrm{N}}\\ & +{\left({\overline{\mathit{B}}}^{\mathrm{N}}\right)}^{\mathrm{T}}\mathit{B}{\overline{\mathit{B}}}^1+{\left({\overline{\mathit{B}}}^1\right)}^{\mathrm{T}}\mathit{B}{\overline{\mathit{B}}}^{\mathrm{N}}+{\left({\overline{\mathit{B}}}^{\mathrm{N}}\right)}^{\mathrm{T}}\mathit{A}{\overline{\mathit{B}}}^1]\mathrm{d}\mathrm{\Omega }.\end{array}

The motion equation of the plates is derived from the nonlinear motion equation of the plate element, as follows:

(47)\mathit{M}\ddot{\mathit{q}}+\left({\overline{\mathit{K}}}_{}^{\mathrm{L}}+{\overline{\mathit{K}}}_{}^{\mathrm{N}}\right)\mathit{q}=\mathit{F}.

These matrices were formed from the element matrices using the finite element algorithm [59].

If structural damping and force vector \mathit{F}=\mathit{F}(t) are included, then the motion equation of the plates (shown in Eq. (47)) can be written in the following form:

(48)\mathit{M}\ddot{\mathit{q}}+\mathit{C}\dot{\mathit{q}}+\left({\overline{\mathit{K}}}_{}^{\mathrm{L}}+{\overline{\mathit{K}}}_{}^{\mathrm{N}}\right)\mathit{q}=\mathit{F}(t),

where \mathit{C}=\alpha \mathit{M}+\beta {\overline{\mathit{K}}}^{\mathrm{L}}. Here, \alpha and \beta are computed as follows [60–64]:

(49)\beta =\frac{2\zeta }{{\omega }_1+{\omega }_2},\alpha =\beta {\omega }_1{\omega }_2,

where \zeta is the drag ratio, and {\omega }_1 and {\omega }_2 are the first two frequencies [59].

The equation above is a second-order nonlinear differential equation. To solve it, we performed Newmark integration combined with the Newton–Raphson algorithm [59], as shown by the flowchart in Fig.3. The solution steps are as follows.

For the constant average acceleration method: \gamma ={\displaystyle \frac{1}{2}};\beta ={\displaystyle \frac{1}{4}}.

Step 1: Initial calculation.

Calculate {\ddot{\mathit{q}}}_0={\displaystyle \frac{{\mathit{R}}_0-\mathit{C}{\dot{\mathit{q}}}_0-{\mathit{f}}_{\mathrm{S}}{\mathit{q}}_0}{\mathit{M}}}, where {\mathit{R}}_0,{\dot{\mathit{q}}}_0, and {\mathit{q}}_0 are the node force vector, initial velocity vector, and displacement vector at time t=0, respectively.

Select \mathrm{\Delta }t.

Calculate a={\displaystyle \frac{\gamma }{\beta }}\mathit{C}+{\displaystyle \frac{1}{\beta {\left(\mathrm{\Delta }t\right)}^{}}}\mathit{M};b={\displaystyle \frac{1}{2\beta }}\mathit{M}+\mathrm{\Delta }t\left({\displaystyle \frac{\gamma }{2\beta }}-1\right)\mathit{C}.

Step 2: Calculation at each time step i.

Calculate \mathrm{\Delta }{\stackrel{⌢}{\mathit{R}}}_i=\mathrm{\Delta }{\mathit{R}}_i+a{\dot{\mathit{q}}}_i+b{\ddot{\mathit{q}}}_i.

Determine the tangential stiffness matrix {\mathit{K}}_i.

Calculate {\stackrel{⌢}{\mathit{K}}}_i={\mathit{K}}_i+{\displaystyle \frac{\gamma }{\beta \mathrm{\Delta }t}}\mathit{C}+{\displaystyle \frac{1}{\beta {\left(\mathrm{\Delta }t\right)}^2}}\mathit{M}.

Determine \mathrm{\Delta }{\mathit{q}}_i from {\stackrel{⌢}{\mathit{K}}}_i and \mathrm{\Delta }{\stackrel{⌢}{\mathit{R}}}_i via the Newton–Raphson iterative procedure.

Calculate \mathrm{\Delta }{\dot{\mathit{q}}}_i={\displaystyle \frac{\gamma }{\beta \mathrm{\Delta }t}}\mathrm{\Delta }{\mathit{q}}_i-{\displaystyle \frac{\gamma }{\beta }}{\dot{\mathit{q}}}_i+\mathrm{\Delta }t\left(1-{\displaystyle \frac{\gamma }{2\beta }}\right){\ddot{\mathit{q}}}_i.

Calculate \mathrm{\Delta }{\ddot{\mathit{q}}}_i={\displaystyle \frac{1}{\beta {\left(\mathrm{\Delta }t\right)}^2}}\mathrm{\Delta }{\mathit{q}}_i-{\displaystyle \frac{1}{\beta \mathrm{\Delta }t}}{\dot{\mathit{q}}}_i-{\displaystyle \frac{1}{2\beta }}{\ddot{\mathit{q}}}_i.

Calculate {\mathit{q}}_{i+1}={\mathit{q}}_i+\mathrm{\Delta }{\mathit{q}}_i;{\dot{\mathit{q}}}_{i+1}={\dot{\mathit{q}}}_i+\mathrm{\Delta }{\dot{\mathit{q}}}_i;{\ddot{\mathit{q}}}_{i+1}={\ddot{\mathit{q}}}_i+\mathrm{\Delta }{\ddot{\mathit{q}}}_i.

Step 3: Repeat the calculation for the next time step in Step 2.

The boundary conditions (BCs) introduced are as follows.

The classic (immovable) boundary edges can be written as

Clamped support (C):

u_0=v_0=w_0={\theta }_x={\theta }_y=0 with all edges.

Simply supported (S):

u_0=w_0={\theta }_x=0 with x=0,a.

v_0=w_0={\theta }_y=0 with y=0,b.

The movable boundary edges can be written as

Simply supported (S1):

w_0={\theta }_x=0 with x=0,a.

w_0={\theta }_y=0 with y=0,b.

Clamped support (C1):

u_0=w_0={\theta }_x={\theta }_y=0 with x=0,a.

v_0=w_0={\theta }_x={\theta }_y=0 with y=0,b.

Free support (F): The degrees of freedom (DOFs) are not equal to zero at the boundary edges.

In this study, we only implemented the immovable BCs.

4 Results and discussion

4.1 Verification studies

First, the dynamic response was compared with the exact solution of Reddy [65] as shown in Fig.4. Herein, an SSSS FGM square plate with geometrical parameters a=b=0.2\mathrm{m},h=a/20 under a uniform sudden load q_0=1\mathrm{M}\mathrm{P}\mathrm{a} is considered. The material properties are as follows: ceramic E_{\mathrm{t}}=151\text{GPa,}{\rho }_{\mathrm{t}}=3000\text{kg/}{\text{m}}^{\text{3}}\text{,}v_{\mathrm{t}}=0.3, metal E_{\mathrm{b}}=70\text{GPa,}{\rho }_{\mathrm{b}}=2707\text{kg/}{\text{m}}^{\text{3}}\text{,} and v_{\mathrm{b}}=0.3. The dimensionless deflection and the dimensionless time are performed based on w^{\ast }=w{\displaystyle \frac{E_{\mathrm{b}}h}{q_0{a}^2}} and t^{\ast }=t\sqrt{{\displaystyle \frac{E_{\mathrm{b}}}{q_0{a}^2}}}. The obtained results converged at the mesh size of 12\times 12 and agreed well with the exact solution [65]. Hereinafter, the mesh size of 12\times 12 (144 elements) is used for the examples.

Next, we consider an SSSS FGM square plate with geometrical dimensions a = b = 1 m, h = a/10. In this example, the mechanical properties are set as follows: ceramic (Si₃N₄) E_{\mathrm{t}}=151\text{GPa}, {\rho }_{\mathrm{t}}=3000\text{kg/}{\text{m}}^{\text{3}}, v_{\mathrm{t}}=0.3; metal (SUS304) E_{\mathrm{b}}=70\text{GPa,}{\rho }_{\mathrm{b}}=2707\text{kg/}{\text{m}}^{\text{3}}, and v_{\mathrm{b}}=0.3. Tab.1 lists a comparison of the nonlinear frequency ratio of SSSS FGM plates (power-law index k = 1) with the results reported by Sundararajan et al. [66] using the first-order Q8 element and those by Balamurugan et al. [67] using the first-order Q9 element. The comparison shows excellent agreement (with the error being less than 2%), where w/h is the flexural amplitude ratio.

The first dimensionless frequencies {\omega }^{\ast }=\omega {\displaystyle \frac{a^2}{h}}\sqrt{{\displaystyle \frac{{\rho }_m}{E_m}}} of the FG-CNTRC square plates are provided in Tab.2 and Tab.3. The obtained results matched closely with the results of the higher-order IGA [68]. Notably, using the Q4 element based on FSDT may not be suitable for analyzing thick FG plates. The MISQ20 element is based on the Q4 element and FSDT; however, it uses the mixed integration smoothed technique to ensure accuracy and reliability. Hence, one can conclude that using the current formula guarantees high accuracy in predicting the nonlinear vibration of plates.

4.2 Nonlinear free vibration

In this section, we consider FG-CNTRC plates with the following material parameters: E_{11}^{\mathrm{C}\mathrm{N}\mathrm{T}}=5.6466\mathrm{T}\mathrm{P}\mathrm{a}, E_2^{\mathrm{C}\mathrm{N}\mathrm{T}}=7.0800\mathrm{T}\mathrm{P}\mathrm{a}, G_{12}^{\mathrm{C}\mathrm{N}\mathrm{T}}=1.9445\mathrm{T}\mathrm{P}\mathrm{a}, v_{12}^{\mathrm{C}\mathrm{N}\mathrm{T}}=0.175, {\rho }_{\mathrm{C}\mathrm{N}\mathrm{T}}=1400\mathrm{k}\mathrm{g}/{\text{m}}^{\text{3}}, E_{}^{\mathrm{m}}=2.5\mathrm{G}\mathrm{P}\mathrm{a}, G_{\mathrm{m}}^{}={\displaystyle \frac{E_{\mathrm{m}}}{2(1+v_{\mathrm{m}})}}, v_{\mathrm{m}}^{}=0.34, and {\rho }_{\mathrm{m}}=1150\text{kg/}{\text{m}}^{\text{3}}.

The first six mode shapes of the UD-type CCCC FG-CNTRC square plate are shown in Fig.5. Its geometrical dimensions are a = b = 1 m and h = a/25, and the volume of the CNT, V_{\mathrm{C}\mathrm{N}\mathrm{T}}^{\ast }=0.12.

Tab.4 and Tab.5 list the simultaneous effects of the four types of CNTs and the geometrical parameters for the nonlinear frequency of the plates. In terms of the CNT distribution type, it appeared that FG-O and FG-X resulted in the maximum and minimum nonlinear frequencies, respectively. Moreover, based on the same BCs and material properties, a plate with a smaller area appeared to be stiffer, thus resulting in a higher nonlinear frequency. As expected, the frequency ratio increased with the flexural amplitude ratio (w/h).

The first six mode shapes of the SSSS FG-CNTRC plate with a center hole are shown in Fig.7. The volume fraction of the CNT, V_{\mathrm{C}\mathrm{N}\mathrm{T}}^{\ast }=0.17 for the FG-V-type distribution, and the geometrical dimensions were a = b = 1 m and h = a/15, as shown in Fig.6. The vibration mode shapes of complex models (such as plates with holes) confirmed the superiority of the numerical method over the analytical method. In this study, we meshed the plate with 256 elements to enhance the smoothness of the mode shapes.

4.3 Nonlinear dynamic response

In this section, we consider the plates subjected to a triangular load P=q_{o}F(t) with q_o=2.1\text{MPa}, as follows:

(51)

where t_{\mathrm{h}\mathrm{d}}=0.03\text{s} is the positive phase duration and the total integration time is t_{\mathrm{c}\mathrm{a}\mathrm{l}}=0.1\text{s}. Notably, a damping ratio \zeta =0.01 was applied to all of the remaining examples.

First, Fig.8 shows the effect of four types of CNT distributions on the nonlinear dynamic response of SSSS FG-CNTRC square plates with geometrical dimensions a=b=0.5\mathrm{m} and h=a/10. Specifically, Fig.8(a) and Fig.8(b) present the nonlinear deflection (w) and velocity (v) response over time, respectively. This proves that the stiffness of the FG-O and FG-X-type plates is the highest and lowest, respectively. Hence, one can conclude that the plate gradually reduces the amplitude of the oscillations if structural damping is considered.

Second, Fig.9 shows the effect of the volume fractions of the CNTs (V_{\mathrm{C}\mathrm{N}\mathrm{T}}^{\ast }) on the dynamic response of the CCCC FG-CNTRC square plates with geometrical dimensions a=b=0.5\text{m} and h=a/15. In this example, the volume of the CNTs, V_{\mathrm{C}\mathrm{N}\mathrm{T}}^{\ast }=0.12,0.17, and 0.28 for the FG-O type distribution. Based on these figures, V_{\mathrm{C}\mathrm{N}\mathrm{T}}^{\ast } increases with the plate stiffness; therefore, the nonlinear deflection and velocity response over time are reduced.

Next, the plot presented in Fig.10 shows the effect of the length-to-thickness ratio (a/h) on the dynamic response of the SCSC FG-CNTRC square plates with geometrical dimensions a=b=0.5\text{m}. The a/h ratios investigated were 10, 25, 50, and 75, and the material parameter V_{\mathrm{C}\mathrm{N}\mathrm{T}}^{\ast }=0.28 for the FG-V type. The result shown in this figure matches our expectation. A plate with greater thickness will become stiffer; therefore, the oscillation amplitude of the thicker plates is smaller in terms of the nonlinear time-history deflection responses.

Finally, Fig.11 shows the effect of the length-to-width ratio (a/b) on the dynamic response of the CFCF FG-CNTRC square plates with geometrical dimensions a=0.5\text{m} (a is fixed) and h=a/20. The a/b ratios investigated were 2, 1.5, 1, and 0.5, and the material parameter V_{\mathrm{C}\mathrm{N}\mathrm{T}}^{\ast }=0.12 for the FG-X type. Based on this figure, one can conclude that the smaller-area plate will become stiffer and thus yield nonlinear deflection and velocity responses over time. In this example, the total investigation time was t_{\mathrm{c}\mathrm{a}\mathrm{l}}=0.3\text{s}.

5 Conclusions

By employing the MISQ20 element, we successfully established a nonlinear motion equation for FG-CNTRC plates. The calculation program was developed using MATLAB software, and based on the coded program, the authors examined the effects of the input parameters on the nonlinear vibration of the plates. The findings of this study are as follows.

1) Using the MISQ20 element provided a higher convergence rate than the Q4 element, which was equal to that of the Q8 element. In addition, the computing time of the MISQ20 element (integrated on the line) was less than that of Q8 because it contained only 20 DOFs, as compared with 40 DOFs in the Q8 element (integrating on the surface). To solve nonlinear problems, a rapid looping process with a high convergence rate is critical.

2) Using the MISQ20 element mitigated the “shear locking” phenomenon without requiring any correction factor. Therefore, employing the MISQ20 element in this study successfully demonstrated its advantages.

3) Based on the same geometric parameters, the FG-O and FG-X-type plates demonstrated the highest and lowest hardness, respectively. Increasing the a/b ratio reduced the plate area, which consequently increased the plate stiffness. To clarify, when the structure size decreases, under the same BC and parameters, the structure will become stronger.

4) The numerical results presented would benefit the calculation and design of FG-CNTRC plates in engineering practice.

5) The proposed program is expected to be further developed for the analysis of nonlinear dynamic problems involving complex structures, which are difficult to solve using analytical methods.

References

[1] Thostenson E T, Ren Z, Chou T W. Advances in the science and technology of carbon nanotubes and their composites: A review. Composites Science and Technology, 2001, 61(13): 1899–1912

[2] Fiedler B, Gojny F H, Wichmann M H G, Nolte M C M, Schulte K. Fundamental aspects of nano-reinforced composites. Composites Science and Technology, 2006, 66(3): 115–125

[3] Cooper C A, Cohen S R, Barber A H, Wagner H D. Detachment of nanotubes from a polymer matrix. Applied Physics Letters, 2002, 81(20): 3873–3875

[4] Barber A H, Cohen S R, Wagner H D. Measurement of carbon nanotube–polymer interfacial strength. Applied Physics Letters, 2003, 82(4): 140–142

[5] GouJMinaieBWangBLiangZZhangC. Computational and experimental study of interfacial bonding of single-walled nanotube reinforced composites. Computational Materials Science, 2004, 31(3−4): 225−236

[6] Frankland S J V, Caglar A, Brenner D W, Griebel M. Molecular simulation of the influence of chemical cross-links on the shear strength of carbon nanotube polymer interfaces. Journal of Physical Chemistry B, 2002, 106(12): 3046–3048

[7] Ma P C, Mo S Y, Tang B Z, Kim J K. Dispersion, interfacial interaction and reagglomeration of functionalized carbon nanotubes in epoxy composites. Carbon, 2010, 48(6): 1824–1834

[8] Coleman J N, Khan U, Blau W J, Gun’ko Y K. Small but strong: A review of the mechanical properties of carbon nanotube–polymer composites. Carbon, 2006, 44(9): 1624–1652

[9] Wagner H D, Lourie O, Feldman Y, Tenne R. Stress-induced fragmentation of multiwall carbon nanotubes in a polymer matrix. Applied Physics Letters, 1998, 72(2): 188–190

[10] Qian D, Dickey E C, Andrews R, Rantell T. Load transfer and deformation mechanisms in carbon nanotube–polystyrene composites. Applied Physics Letters, 2000, 76(20): 2868–2870

[11] Odegard G M, Gates T S, Nicholson L M, Wise K E. Equivalent-continuum modeling of nano-structured materials. Composites Science and Technology, 2002, 62(14): 1869–1880

[12] Odegard G M, Gates T S, Wise K E, Park C, Siochi E J. Constitutive modeling of nanotube-reinforced polymer composites. Composites Science and Technology, 2003, 63(11): 1671–1687

[13] Liu Y J, Chen X L. Evaluations of the effective material properties of carbon nanotube-based composites using a nanoscale representative volume element. Mechanics of Materials, 2003, 35(1−2): 69–81

[14] Hu N, Fukunaga H, Lu C, Kameyama M, Yan B. Prediction of elastic properties of carbon nanotube reinforced composites. Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences, 2005, 461(2058): 1685–1710

[15] Frankland S J V, Harik V M, Odegard G M, Brenner D W, Gates T S. The stress–strain behaviour of polymer–nanotube composites from molecular dynamics simulation. Composites Science and Technology, 2003, 63(11): 1655–1661

[16] Griebel M, Hamaekers J. Molecular dynamics simulations of the elastic moduli of polymer–carbon nanotube composites. Computer Methods in Applied Mechanics and Engineering, 2004, 193(17−20): 1773–1788

[17] Wernik J M, Meguid S A. Multiscalemodeling of the nonlinear response of nanoreinforced polymers. Acta Mechanica, 2011, 217(1−2): 1–16

[18] Wuite J, Adali S. Deflection and stress behaviour of nanocomposite reinforced beams using a multiscale analysis. Composite Structures, 2005, 71(3−4): 388–396

[19] Vodenitcharova T, Zhang L C. Bending and local buckling of a nanocomposite beam reinforced by a single-walled carbon nanotube. International Journal of Solids and Structures, 2006, 43(10): 3006–3024

[20] Ray M C, Batra R C. A single-walled carbon nanotube reinforced 1–3 piezoelectric composite for active control of smart structures. Smart Materials and Structures, 2007, 16(5): 1936–1947

[21] Formica G, Lacarbonara W, Alessi R. Vibrations of carbon nanotube-reinforced composites. Journal of Sound and Vibration, 2010, 329(10): 1875–1889

[22] Arani A, Maghamikia S, Mohammadimehr M, Arefmanesh A. Buckling analysis of laminated composite rectangular plates reinforced by SWCNTS using analytical and finite element methods. Journal of Mechanical Science and Technology, 2011, 25(3): 809–820

[23] Shen S H. Nonlinear bending of functionally graded carbon nanotube-reinforced composite plates in thermal environments. Composite Structures, 2009, 91(1): 9–19

[24] Wang Z X, Shen H S. Nonlinear vibration of nanotube-reinforced composite plates in thermal environments. Computational Materials Science, 2011, 50(8): 2319–2330

[25] Wang Z X, Shen H S. Nonlinear vibration and bending of sandwich plates with nanotube-reinforced composite face sheets. Composites. Part B, Engineering, 2011, 43(2): 411–421

[26] Zhang L, Cui W, Liew K. Vibration analysis of functionally graded carbon nanotube reinforced composite thick plates with elastically restrained edges. International Journal of Mechanical Sciences, 2015, 103: 9–21

[27] Natarajan S, Haboussi M, Manickam G. Application of higher-order structural theory to bending and free vibration analysis of sandwich plates with CNT reinforced composite facesheets. Composite Structures, 2014, 113: 197–207

[28] Sankar A, Natarajan S, Haboussi M, Ramajeyathilagam K, Ganapathi M. Panel flutter characteristics of sandwich plates with CNT reinforced facesheets using an accurate higher-order theory. Journal of Fluids and Structures, 2014, 50: 376–391

[29] Sankar A, Natarajan S, Zineb T B, Ganapathi M. Investigation of supersonic flutter of thick doubly curved sandwich panels with CNT reinforced facesheets using higher-order structural theory. Composite Structures, 2015, 127: 340–355

[30] Sankar A, Natarajan S, Merzouki T, Ganapathi M. Nonlinear dynamic thermal buckling of sandwich spherical and conical shells with CNT reinforced facesheets. International Journal of Structural Stability and Dynamics, 2017, 17(9): 1750100

[31] Rodrigues J, Natarajan S, Ferreira A, Carrera E, Cinefra M, Bordas S. Analysis of composite plates through cell-based smoothed finite element and 4-noded mixed interpolation of tensorial components techniques. Computers & Structures, 2014, 135: 83–87

[32] HughesT J RCottrellJ ABazilevsY. Isogeometric analysis: CAD, finite elements, NURBS, exact geometry and mesh refinement. Computer Methods in Applied Mechanics and Engineering, 2005, 194(39−41): 4135−4195

[33] BordenM JScottM AEvamsJ AHughesT J R. Isogeometric finite element data structures based on Bezier extraction of NURBS. International Journal for Numerical Methods in Engineering, 2011, 87 (1−5): 15−47

[34] Phung-Van P, Thai C H, Nguyen-Xuan H, Abdel-Wahab M. An isogeometric approach of static and free vibration analyses for porous FG nanoplates. European Journal of Mechanics. A, Solids, 2019, 78: 103851

[35] Thanh C L, Nguyen T N, Vu T H, Khatir S, Abdel Wahab M. A geometrically nonlinear size-dependent hypothesis for porous functionally graded micro-plate. Engineering with Computers, 2022, 38(Suppl 1): 449–460

[36] Cuong-Le T, Nguyen K D, Hoang-Le M, Sang-To T, Phan-Vu P, Wahab M A. Nonlocal strain gradient IGA numerical solution for static bending, free vibration and buckling of sigmoid FG sandwich nanoplate. Physica B, Condensed Matter, 2022, 631: 413726

[37] PhamQ HNguyenP CTranV KLieuQ XTranT T. Modified nonlocal couple stress isogeometric approach for bending and free vibration analysis of functionally graded nanoplates. Engineering with Computers, 2022: 1−26

[38] Pham Q H, Nguyen P C, Tran T T. Dynamic response of porous functionally graded sandwich nanoplates using nonlocal higher-order isogeometric analysis. Composite Structures, 2022, 290: 115565

[39] Guo H, Zheng H. The linear analysis of thin shell problems using the numerical manifold method. Thin-walled Structures, 2018, 124: 366–383

[40] Guo H, Zheng H, Zhuang X. Numerical manifold method for vibration analysis of Kirchhoff’s plates of arbitrary geometry. Applied Mathematical Modelling, 2019, 66: 695–727

[41] Zhuang X, Guo H, Alajlan N, Zhu H, Rabczuk T. Deep autoencoder based energy method for the bending, vibration, and buckling analysis of Kirchhoff plates with transfer learning. European Journal of Mechanics. A, Solids, 2021, 87: 104225

[42] Guo H, Zhuang X, Rabczuk T. A deep collocation method for the bending analysis of Kirchhoff plate. Computers Materials & Continua, 2021, 59(2): 433–456

[43] Samaniego E, Anitescu C, Goswami S, Nguyen-Thanh V M, Guo H, Hamdia K, Zhuang X, Rabczuk T. An energy approach to the solution of partial differential equations in computational mechanics via machine learning: Concepts, implementation and applications. Computer Methods in Applied Mechanics and Engineering, 2020, 362: 112790

[44] Ganapathi M, Varadan T K, Sarma B S. Nonlinear flexural vibrations of laminated orthotropic plates. Computers & Structures, 1991, 39(6): 685–688

[45] Kant T, Kommineni J R. Large amplitude free vibration analysis of cross-ply composite and sandwich laminates with a refined theory and C⁰ finite elements. Computers & Structures, 1994, 50(1): 123–134

[46] AnilK D. Large amplitude free vibration analysis of composite plates by finite element method. Thesis for the Master’s Degree. Rourkela: National Institute of Technology, 2010

[47] JavedA. Dynamic stability of delaminated cross ply composite plates and shells. International Journal of Mechanical Sciences, 1998, 40(8): 805–823

[48] Parhi P K, Bhattacharyya S K, Sinha P K. Hygrothermal effects on the dynamic behaviour of multiple delamated composite plates and shells. Journal of Sound and Vibration, 2001, 248(2): 195–214

[49] Providas E, Kattis M A. An assessment of two fundamental flat triangular shell elements with drilling rotations. Computers & Structures, 2000, 77(2): 129–139

[50] Groenwold A A, Slander N. An efficient 4-node 24 DOF thick shell finite element with 5-point quadrature. Engineering Computations, 1995, 12(8): 723–747

[51] Choi C K, Lee T Y. Efficient remedy for membrane locking of 4-node flat shell elements by non-conforming modes. Computer Methods in Applied Mechanics and Engineering, 2003, 192(16–18): 1961–1971

[52] Pimpinelli G. An assumed strain quadrilateral element with drilling degrees of freedom. Finite Elements in Analysis and Design, 2004, 41(3): 267–283

[53] Thai-Hoang C, Nguyen-Thanh N, Nguyen-Xuan H, Rabczuk T. An alternative alpha finite element method with discrete shear gap technique for analysis of laminated composite plates. Applied Mathematics and Computation, 2011, 217(17): 7324–7348

[54] Phan-Dao H, Nguyen-Xuan H, Thai-Hoang C, Nguyen-Thoi T, Rabczuk T. An edge-based smoothed finite element method for analysis of laminated composite plates. International Journal of Computational Methods, 2013, 10(1): 1340005

[55] Nguyen-Thanh N, Rabczuk T, Nguyen-Xuan H, Bordas S. A smoothed finite element method for shell analysis. Computer Methods in Applied Mechanics and Engineering, 2008, 198(2): 165–177

[56] Nguyen-XuanHRabczukTBordasSDebongnieJ F. A smoothed finite element method for plate analysis. Computer Methods in Applied Mechanics and Engineering, 2008, 197(13−16): 1184−1203

[57] Shen H S. Postbuckling of nanotube-reinforced composite cylindrical shells in thermal environments, Part I: Axially loaded shells. Composite Structures, 2011, 93(8): 2096–2108

[58] Nguyen-Van H, Nguyen-Hoai N, Chau-Dinh T, Nguyen-Thoi T. Geometrically nonlinear analysis of composite plates and shells via a quadrilateral element with good coarse-mesh accuracy. Composite Structures, 2014, 112: 327–338

[59] ReddyJ N. Mechanics of Laminated Composite Plates and Shells: Theory and Analysis. Boca Raton: CRC Press, 2004

[60] ChopraA K. Dynamics of Structures: Theory and Applications to Earthquake Engineering. Upper Saddle River: Pearson Prentice Hall, 2007

[61] Nguyen P C, Kim S E. Investigating effects of various base restraints on the nonlinear inelastic static and seismic responses of steel frames. International Journal of Non-linear Mechanics, 2017, 89: 151–167

[62] TranTTPhamQHNguyen-ThoiT. Dynamic analysis of sandwich auxetic honeycomb plates subjected to moving oscillator load on elastic foundation. Advances in Materials Science and Engineering, 2020, 6309130

[63] Nguyen H N, Canh T N, Thanh T T, Ke T V, Phan V D, Thom D V. Finite element modelling of a composite shell with shear connectors. Symmetry, 2019, 11(4): 527

[64] Tran T T, Tran V K, Le P B, Phung V M, Do V T, Nguyen H N. Forced vibration analysis of laminated composite shells reinforced with graphene nanoplatelets using finite element method. Advances in Civil Engineering, 2020, 2020: 1471037

[65] Reddy J. Analysis of functionally graded plates. International Journal for Numerical Methods in Engineering, 2000, 47(1–3): 663–684

[66] Sundararajan N, Prakash T, Ganapathi M. Nonlinear free exural vibrations of functionally graded rectangular and skew plates under thermal environments. Finite Elements in Analysis and Design, 2005, 42(2): 152–168

[67] Balamurugan V, Ganapathi M, Varadan T. Nonlinear dynamic instability of laminated composite plates using finite element method. Computers & Structures, 1996, 60(1): 125–130

[68] Phung-Van P, Abdel-Wahab M, Liew K, Bordas S, Nguyen-Xuan H. Isogeometric analysis of functionally graded carbon nanotube-reinforced composite plates using higher-order shear deformation theory. Composite Structures, 2015, 123: 137–149

RIGHTS & PERMISSIONS

Higher Education Press