1 Introduction

In recent years, the mathematical modeling and numerical simulation of small-scale structures, especially carbon nanotubes (CNTs), have become noteworthy new subjects in materials and engineering science. Hence, many investigations have been performed to study various aspects of these structures^[1-4].

Recently, a new class of advanced composites, namely, functionally graded carbon nanotube-reinforced composite (FG-CNTRC), has received a wide range of attention in industrial, engineering, and research societies due to excellent and exclusive mechanical properties as compared with those of conventional composites. Owing to their excellent mechanical properties^[5-6], FG-CNTRCs can be employed in several applications such as reinforcing composites, electronic devices, new nano- and micro-electromechanical systems, and many others. Hence, the mathematical modeling of beam-, plate-, and shell-like structures made of FG-CNTRCs has drawn a lot of attention from researchers to understand their linear and nonlinear static and dynamic mechanical behaviors^[7-16].

1.1 Literature review on the analysis of the mechanical behavior of perfect FG-CNTRC structures

For FG-CNTRC beams, the geometrically nonlinear forced vibration of FG-CNTRC Timoshenko beams with various boundary conditions was numerically examined by Ansari et al.^[9]. Lin and Xiang^[17] performed a numerical analysis on the nonlinear free vibration of first-order and third-order shear deformable FG-CNTRC beams. Using a numerical solution approach, the free vibration and axial buckling of sandwich Timoshenko beams reinforced with nanocomposite face sheets were analyzed by Wu et al.^[18]. The theory of elasticity and the state-space differential quadrature method were utilized by Alibeigloo and Liew^[19] to simulate the free vibration and bending characteristics of piezoelectric FG-CNTRC beams. For FG-CNTRC plates and in the category of linear study, Zhang et al.^[20] utilized the element-free approach to analyze the buckling of FG-CNTRC skew plates. Moreover, with the aid of the state-space Levy method, the free vibration of third-order shear deformable FG-CNTRC plates was examined by Zhang et al.^[21]. Zhu et al.^[22] attempted to study the bending and free vibration of first-order shear deformable FG-CNTRC rectangular plates using the finite element approach. In the category of nonlinear analysis, the nonlinear bending of FG-CNTRC skew, rectangular, and straight-sided quadrilateral plates embedded in the Pasternak elastic foundation was numerically investigated by utilizing the element-free improved moving least-squares (IMLS) Ritz method^[23-25]. Using a numerical solution strategy, Ansari et al.^[26] tried to examine the nonlinear forced vibration of FG-CNTRC rectangular plates subject to a harmonic excitation transverse load. Also, to accurately simulate the nonlinear primary resonance of thick and moderately thick FG-CNTRC plates with various boundary conditions without using the shear correction factor, Ansari and Gholami^[27] developed a nonlinear third-order shear deformable plate model by making use of the von Kármán assumptions and Reddy's third-order shear deformation plate theory. Various types of the periodic and chaotic motions of FG-CNTRC plate subject to the combination of parametric and forcing loadings were analyzed by Guo and Zhang^[28]. Wang and Shen^[29] carried out nonlinear vibration of embedded FG-CNTRC plates in the thermal environment using an analytical solution method. For FG-CNTRC shells, an analytical solution procedure was provided by Ansari et al.^[30] to obtain the nonlinear secondary equilibrium postbuckling path of FG-CNTRC cylindrical shells with piezoelectric layers. Using a numerical solution procedure, Lei et al.^[31] analyzed the dynamic stability of FG-CNTRC cylindrical panels under the periodic compressive axial forces. In the context of higher-order shear deformation shell theory and applying the von Kármán-type of kinematic nonlinearity, an analytical solution was utilized by Shen^[32] to examine the nonlinear thermal postbuckling of FG-CNTRC cylindrical shells. Recently, Ansari and Torabi^[33] investigated the free vibration and axial buckling of FG-CNTRC conical shells.

1.2 Literature review on the analysis of the initially imperfect FG-CNTRC structures

It must be remarked that the previous investigations were restricted to static and dynamic characteristics of perfect straight FG-CNTRC structures only. Initial geometric imperfections may be created in the structures due to the inappropriate manufacturing process or during their service life. Furthermore, many structures with initial geometric imperfections can be found that are utilized in various industry and engineering applications. The initial geometric imperfections may result in considerable effects on the mechanical behaviors of structures. In the real conditions, the initial geometric imperfections may arise in different shapes. However, in many studies, the initial imperfection in the structures was supposed as the mode shapes of the buckling or free vibration of system. Studies on the impact of geometric imperfection on the mechanical behaviors of FG-CNTRC systems are still limited in number. In this regard, by employing the Timoshenko beam theory and Ritz approach, Wu et al.^[34] tried to analyze the nonlinear free vibration of imperfect FG-CNTRC beams. They investigated the effect of the geometric imperfections with different shapes including sine type, global, and localized imperfections on the nonlinear frequency-amplitude curve of FG-CNTRC beams. Also, the imperfection sensitivity of mechanical and thermal postbuckling of FG-CNTRC beams was examined by Wu et al.^[35-36]. Thang et al.^[37] analytically studied the nonlinear postbuckling of fully simply-supported imperfect FG-CNTRC plates based on the Kirchhoff plate theory.

1.3 Contributions of this study to the field

To the authors' best knowledge, the geometrically nonlinear resonant dynamics of initially imperfect FG-CNTRC rectangular plates is still unexplored. Moreover, in most of the above-mentioned studies, the Kirchhoff and first-order shear deformation plate theories have been used which cannot correctly predict the mechanics of thick and moderately thick plates. Therefore, this study, for the first time, aims to develop the mathematical formulation of initially imperfect shear deformable FG-CNTRC rectangular plates. The proposed model includes the shear deformation effect and rotary inertia without the need for any shear correction factor. Then, on the basis of the proposed model, this work examines the impact of initial geometric imperfection on the geometrically nonlinear dynamical behaviors of shear deformable FG-CNTRC rectangular plates with various edge conditions under harmonic excitation transverse load. The initial imperfection is assumed to be as the first mode of free vibration of perfect FG-CNTRC plate corresponding to each edge support and given initial deflection. The rule of mixture is utilized to calculate the effective material properties of FG-CNTRC rectangular plates. Within the framework of the parabolic shear deformation plate theory and applying the von Kármán hypotheses, Hamilton's principle is utilized to achieve the geometrically nonlinear mathematical formulation of initially imperfect FG-CNTRC plates. The developed plate model includes the effect of initial geometric imperfection, rotary inertia, and transverse shear deformation without using the shear correction factor. These advantages enable us to analyze the nonlinear dynamical characteristics of initially imperfect thick and moderately thick FG-CNTRC rectangular plates using the developed shear deformable plate model. Afterwards, to analyze the imperfection sensitivity of the nonlinear resonant dynamics of imperfect FG-CNTRC plates, a multistep solution methodology including the generalized differential quadrature (GDQ) method, the Galerkin approach, and the time periodic discretization (TPD) scheme is employed to convert the geometrically nonlinear partial differential equations into a set of parameterized nonlinear equations. Then, a combination of the pseudo-arc length iterative technique and the modified Newton-Raphson method is adopted to solve the set of nonlinear algebraic equations to track the frequency- and force-response curves of initially imperfect FG-CNTRC rectangular plates with various edge conditions. Several numerical results are presented to understand the effects of initial geometric imperfection, geometrical parameters, and edge conditions on nonlinear resonant dynamical responses of the imperfect FG-CNTRC rectangular plates.

2 Model development and solution procedure

The nonlinear mathematical formulations of equations of initially imperfect FG-CNTRC rectangular plates are provided in this section. In this regard, the parabolic shear deformation plate theory in conjugation with the von Kármán assumptions is utilized to extract the nonlinear equations of motion of imperfect FG-CNTRC rectangular plates using Hamilton's principle.

An FG-CNTRC rectangular plate with length a, width b, and thickness h subject to an initial geometric imperfection is considered. The plate is composed of the mixture of a matrix and single-walled carbon nanotubes (SWCNTs). It is assumed that the initial imperfection w^*(x₁, x₂) is in the form of the first vibration mode shape of perfect FG-CNTRC rectangular plate corresponding to each edge support with a given value of the deflection γ. Also, the FG-CNTRC rectangular plate is excited externally by a uniformly-distributed harmonic excitation transverse load per unit length , where f₀, , and t indicate the forcing amplitude, the excitation frequency, and the time, respectively. Introducing a Cartesian coordinate system with the origin at one corner (0≤ x₁ ≤ a, 0≤ x₂ ≤ b, -h/2≤ x₃ ≤ h/2), the in-plane displacements of each point of middle-plane along the axes x_I (I=1, 2) are denoted by as u_I. Moreover, the transverse displacement of an arbitrary point on the mid-plane is symbolized by w. Also, the rotations of the middle surface normal at z=0 about the x₂- and x₁-axes are denoted by ψ₁ and ψ₂, respectively.

2.1 Displacement vector, nonlinear strain-displacement relations, and constitutive equations

Utilizing the parabolic shear deformation plate theory^[38-39], the displacement vector u of FG-CNTRC plate can be expressed as

(1)

where is the shape function of shear deformation. It is not unique, and many choices are available. Generally, other forms of shape functions can be selected so that zero traction boundary conditions on the lower and upper surfaces of plates are satisfied. Herein, it is assumed that ^[40]. Furthermore, it should be noted that compared with the classical Kirchhoff plate theory and Mindlin's first-order shear deformation plate theory, although higher-order shear deformation plate theories involve more algebraic complexity and computational efforts, they do not need any shear correction factor, and hence, one can better represent the kinematics and predict stress distributions of structures, especially for thick and moderately thick ones^[41]. Also, it is noted that some complexities of higher-order shear deformation theories can be reduced using the two- and four-variable refined shear deformation theories. Recently, the higher-order shear deformation theories including the third-order, parabolic, trigonometric, hyperbolic, exponential, and four-variable refined shear deformation theories have been extensively utilized in mathematical modeling of mechanical structures^[42-46].

Moreover, ũ_i identifies the elements of displacement of a given point of FG-CNTRC plate along the axes x_i. It should be noted that the symbol comma in the mathematical relations indicates the partial differentiation with respect to the space coordinates. Furthermore, the capital Latin indices I and J can possess the given values, while i, j=1, 2, 3 and α, β =1, 2.

According to the von Kármán assumptions and applying the initial geometric imperfection, the geometrically nonlinear strain-displacement relations on the basis of the displacement vector defined in Eq. (1) are given by

(2)

where

(3)

Now, on the basis of the linear elasticity, the stress-strain relations for FG-CNTRC rectangular plates are expressed as

(4)

in which the elastic moduli Q_IJ (I, J=1, 2, 4, 5, 6) of the FG-CNTRCs are calculated as follows:

(5)

where E_IJ, G_IJ, and ν_IJ are the effective Young's moduli, shear moduli, and Poisson's ratios of CNTRCs, respectively, and can be approximately calculated via the rule of mixture^[47-48]. By making use of this rule, one can write the following relations^[49]:

(6a)

(6b)

(6c)

(6d)

(6e)

in which V_cnt denotes the volume fraction of CNTs. As mentioned before, the considered rectangular plate consists of an isotropic matrix and SWCNTs. The CNTs can be uniformly distributed (UD) of functionally graded through the thickness of the plate. The functionally graded distribution can be considered as three types, namely, FGX, FGO, and FGA. For each kind of considered CNT distribution, the volume fraction of CNTs can be computed as follows:

(7a)

(7b)

(7c)

in which

(8)

In Eqs. (6) and (8), the superscripts and subscripts m and cnt, respectively, indicate the matrix and CNT phases. The parameters Λ_cnt, ρ^cnt, and ρ^m are the mass fraction of CNT, the CNT density, and the density of matrix, respectively. Moreover, E, G, and ν represent Young's modulus, the shear modulus, and Poisson's ratio, respectively. η_j (j=1, 2, 3) stand for the CNT efficiency parameters, indicating the load transfer between the CNTs and polymeric parts and any small scale effect. Also, V_m is the volume fraction of matrix, which can be calculated as V_m =1-V_cnt.

2.2 Elastic strain energy, kinetic energy, and external work

On the basis of the Hamilton's principle, one can express the following relation:

(9)

in which δΠ_T is the variation of kinetic energy, and δΠ_s indicates the variation of elastic energy. Also, δΠ_P is the variation of external work due to q(t) in the out-of-plane direction.

Using the components of strain tensor defined in Eq. (2), the variation of elastic strain energy of FG-CNTRC rectangular plate can be expressed as

(10)

where N_αβ, M_αβ, P_αβ, and Q_α stand for the in-plane force resultants, bending moment resultants, higher-order bending moments, and transverse forces, respectively. The aforementioned quantities due to the stress elements σ_αβ and σ_α3 are calculated as

(11)

Making use of Eqs. (2) and (4) and taking the integration of the thickness, the resultants defined in Eq. (11) can be calculated as

(12a)

(12b)

(12c)

(12d)

in which γ₁₂ⁱ =2ε₁₂ⁱ (i=0, 1, 2) and the stiffness coefficients appeared in Eq. (12) are defined as

(13)

Now, substituting the relation (3) into Eq. (11) and subsequently applying the calculus of variations, the variation of elastic strain energy of imperfect FG-CNTRC rectangular plate is rewritten as follows:

(14)

where

(15)

Also, n_α indicate the direction cosines of the outward unit normal to the boundary of the middle plane.

Furthermore, using Eq. (1), the kinetic energy of imperfect FG-CNTRC rectangular plate is expressed as

(16)

where the differentiation with respect to the time is signified by dot. Also, I_J (J=0, 1, ..., 5) are calculated as

(17)

Using the calculus of variation, by taking integration over the time and performing the integration by parts, the variation of kinetic energy of system defined in Eq. (16) can be expressed as follows:

(18)

Moreover, the variation of external work due to the external uniformly-distributed harmonic excitation load per unit area applying on the imperfect FG-CNTRC rectangular system can be written as

(19)

2.3 Nonlinear equations of motion of imperfect FG-CNTRC plates

By inserting Eqs. (14), (18), and (19) into Eq. (9) and applying the lemma of calculus of variations, the nonlinear mathematical formulations related to the motion of initially imperfect FG-CNTRC rectangular plates on the basis of parabolic shear deformation plate theory can be achieved as

(20a)

(20b)

(20c)

and the mathematical expressions for essential and natural boundary conditions are achieved as

(21a)

(21b)

(21c)

(21d)

(21e)

(21f)

By choosing x₁ =x, x₂ =y, u₁ =u, u₂ =v, ψ₁ =ψ_x, and ψ₂ =ψ_y and substituting the relations defined in Eq. (12) into Eq. (20), the nonlinear equations of motions in terms of elements of displacements can be written as follows:

(22a)

(22b)

(22c)

(22d)

(22e)

The nonlinear parameters Z_i (i=1, 2, ..., 5) and linear components in the above equations are defined in Appendix A.

In this study, the nonlinear resonance of imperfect FG-CNTRC plates with three types of boundary conditions, namely, fully clamped (CCCC), fully simply supported (SSSS), and two opposite edges clamped and remaining edges simply-supported (CSCS) are analyzed. Utilizing the achieved relations in Eq. (21), these edge conditions are mathematically stated as follows:

For SSSS,

(23)

For CCCC,

(24)

For CSCS,

(25)

3 Method of solution

To predict the nonlinear resonance dynamics of initially imperfect FG-CNTRC rectangular plates, an efficient multistep numerical solution method is used. This scheme includes the GDQ method, the numerical-based Galerkin approach, and the TPD scheme as well as the combination of the pseudo-arc length iterative technique and the modified Newton-Raphson method. Since the utilized solution methodology has been frequently used by the same authors in previously published works^[50-52], for the sake of brevity, the details of the utilized procedure as well as the discretized counterparts of nonlinear equations are not provided herein.

Briefly, the outline of the used solution procedure is as follows:

(ⅰ) The GDQ approach^[53] is implemented to discretize the nonlinear equations of motion and corresponding end conditions. Applying the GDQ method to Eq. (22) results in the following discretized relation:

(26)

in which represents the displacement vector containing 5 nm components of generalized coordinates (n and m denote the number of mesh grid points along the x- and y-directions, respectively). Also, K, M, F, and K_nl(X) stand for the stiffness matrix, mass matrix, force vector, and nonlinear stiffness vector, respectively.

(ⅱ) By neglecting the nonlinear terms and initial imperfection and considering the solution of free vibration problem as (ω_L denotes the linear frequency), the eigenvalue problem corresponding to the linear free vibration problem is achieved as follows:

(27)

which is solved easily to obtain the mode shapes corresponding to different edge conditions. For each boundary condition, the initial imperfection w^* (x₁, x₂) of the FG-CNTRC rectangular plate is assumed as the form of the first vibration mode shape of perfect plate with specified value of deflection γ.

(ⅲ) Substituting the defined initial imperfection into the discretized formulation and then neglecting the nonlinear terms, the fundamental frequencies and associated mode shapes of initially imperfect FG-CNTRC plates can be achieved via solving the corresponding eigenvalue problem. By choosing the first k mode shapes, one can express the linear response of imperfect FG-CNTRC rectangular as follows:

(28)

in which Φ and q are the Galerkin's basis function and the reduced generalized coordinate vector, respectively. Substituting the relation (28) into Eq. (26) results in the following residual:

(29)

(ⅳ) Making use of the obtained mode shapes and introducing the matrix operator G=Φ^TS (S is the integral operator), the numerical Galerkin method is employed to decrease the dimension of the discretized nonlinear equations and subsequently obtain a Duffing-type set of time-dependent ordinary differential equations of motion as follows:

(30)

in which

(31)

(ⅴ) Assume that the non-conservative viscous damping forces are ^[54] (c is the damping coefficient). The Duffing-type equations in the time domain can be discretized by adopting the TPD scheme to transform them to a set of nonlinear algebraic parameterized equations. Therefore, one can obtain the following relation:

(32)

where and the even number N_t signify the time differentiation operators and the discrete grid points in the time domain, respectively. Furthermore, Q stands for the nodal values of q(τ^*).

At last, one can use a vectorized form φ_vec to write Eq. (32) in the following form:

(33)

(ⅵ) A combination of the pseudo-arc length iterative technique and the modified Newton-Raphson method^[55] is employed to solve the set of nonlinear algebraic equations and to track the frequency-amplitude and forcing-amplitude curves of initially imperfect FG-CNTRC rectangular plates.

4 Results and discussion

Numerical results in the form of graphs are provided in this section for the nonlinear resonant dynamics of imperfect FG-CNTRC plates to clarify the impact of initial geometric imperfection on the nonlinear frequency- and force-responses curves. In this regard, polymethylmethacrylate, referred to as PMMA, with the following material properties at the ambient temperature, is considered as a matrix^[49]

Moreover, for single-walled (10, 10) armchair CNTs, the material properties reported by Shen and Zhang^[56] achieved by MD simulations are utilized in this study,

Also, in order to calculate the CNT efficiency parameters appeared in Eq. (1), Young's moduli (E₁₁, E₂₂) and shear modulus (G₁₂) of CNTRC estimated via the extended rule of mixture must be matched with those obtained using the MD simulations^[57-58]. These CNT efficiency parameters associated with three considered CNT volume fractions are obtained as

It should be noted that for the presented results, the non-dimensional frequencies are denoted by , where A₁₁₀ and I₀₀ denote the values of A₁₁ and I₀ associated with the homogeneous matrix plate. Moreover, the non-dimensional applied forcing, the maximum vibration amplitudes, and the damping coefficient are, respectively, defined as f=f₀a²/hA₁₁₀, w_max → w_max /h, and to describe the numerical results.

First of all, to have a convergence test and verify the accuracy of the presented results, in the special case, the frequency parameters of perfect CCCC and SSSS FGX-CNTRC plates corresponding to various numbers of mesh grid points are provided and compared with those existing in the literature (see Table 1). It is assumed that h/b=0.1, a/b=1, and V_cnt^* =0.17. It can be seen that the converged results are achieved when the FGX-CNTRC plate is discretized with grid points higher than 17× 17. Furthermore, the present results are in good agreement with those provided by Selim et al.^[59].

Figures 1 and 2 depict the differences between the nonlinear frequency-amplitude response curves of fully simply-supported imperfect and perfect FGX-CNTRC plates, respectively, where w_max denotes the non-dimensional maximum amplitude and Ω/ω_L denotes the frequency ratio. On the basis of utilized numerical solution procedure, the stable and unstable solution branches as well as the limit point bifurcations are numerically determined. It can be seen that the nonlinear frequency-response curve of perfect system illustrates a hardening-type nonlinear behavior. Furthermore, in addition to two point bifurcations, it includes only an unstable branch, whereas for the considered system with initial geometric imperfection γ =0.4, the frequency-response system shows a softening-type behavior at small vibration amplitude and then follows by a hardening-type behavior at higher vibration amplitude.

Moreover, it contains four point bifurcations and two unstable branches. The aforementioned differences are due to the degree of importance of quadratic and cubic nonlinear components into formulations. For an imperfect FGX-CNTRC plate, since the quadratic nonlinear terms at small vibration amplitude are dominant, the system displays a softening-type nonlinearity. For a perfect system as well as the imperfect plate at higher vibration amplitude, the cubic nonlinear counterparts are major, resulting in a hardening nonlinear behavior. However, it should be indicated that regarding the amount of initial geometric imperfection, edge supports, and geometric parameters, no softening behavior may be displayed in the frequency-amplitude curve of FG-CNTRC plates. In this regard, the nonlinear frequency-amplitude curves for various length-to-thickness ratios are illustrated in Fig. 3. At small vibration amplitudes, it can be inferred that the increase in the length-to-thickness ratio results in decreasing the softening-type behavior of plate such that for the plate with a high length-to-thickness ratio, no softening behavior can be displayed in the frequency-amplitude curve. Also, the reduction in the length-to-thickness ratio causes to intensifying the softening behavior of system at high vibration amplitudes. Furthermore, for all types of edge conditions, it can be seen that the peak of frequency-response curve increases with an increase in the length-to-thickness ratio of rectangular plates.

Figure 4 shows the force-response curves of fully simply-supported FGX-CNTRC plates with initial geometric imperfection γ =0.4 and Ω /ω_L =0.95, where f denotes the forcing amplitude. It can be found that the curve contains two unstable branches and two jumps as well as four limit point bifurcations, while as indicated by Ansari and Gholami^[27], for Ω /ω_L < 1, no jump and limit bifurcation point can be displayed in the forcing-amplitude curve.

For more details, the impact of length-to-thickness ratio (a/h) on the forcing-amplitude curves is examined in Fig. 5. Corresponding to various boundary conditions and geometric parameters, these unstable solutions are decreased to zero, one, or two branch/branches. For example, in the case of FGX-CNTRC plate with clamped edge condition and a/h=16, the curve illustrates only a stable response without any jump and bifurcation, while for a/h=10, an unstable branch and two limit point bifurcations are displayed with varying the forcing amplitude. However, for a fully simply-supported FGX-CNTRC plate with all considered a/h, the system displays four limit point bifurcations. However, it can be found that in the plate with smaller length-to-thickness ratios, the energy transfer mechanism is stronger.

Figures 6 and 7 investigate the effect of initial geometric imperfection on the nonlinear frequency- and force-responses curves of FGX-CNTRC rectangular plates with various boundary conditions. Moreover, the linear frequency of system associated with each value of initial imperfection is provided for a comparison. It can be found that an increase in the initial imperfection results in increasing the linear natural frequency of FGX-CNTRC rectangular plate. Moreover, as illustrated in the frequency-amplitude response curve, one can find that for perfect plates or FG-CNTRC plates with small initial imperfection, the frequency-amplitude response curve illustrates only hardening-type nonlinearity. However, the hardening-type nonlinearity of system weakens, and the type of nonlinearity in FGX-CNTRC rectangular plates transforms from hardening manner to softening treatment, especially at small vibration amplitude as the initial geometric imperfection increases. Moreover, it is clearly illustrated that compared with the plates with CCCC and CSCS ends, more increase in the softening behavior of the SSSS FGX-CNTRC plate can be found when the initial geometric increases.

Also, for the case of Ω /ω_L < 1, Fig. 7 is the effect of magnitude of initial geometric imperfection on the forcing-amplitude curve of plates with different boundary conditions. This figure shows no limit bifurcation point, unstable branch, and jump phenomena in the nonlinear forcing-amplitude curve of a perfect FGX-CNTRC plate. Hence, raising the forcing amplitude leads to gradually increasing the vibration amplitude of system. However, regarding the type of end conditions, the value of initial imperfection and geometry of plate, either two or four limit point bifurcations and one or two unstable solutions may occur in the nonlinear forcing-amplitude curve.

5 Concluding remarks

The aim of this study is to examine the impact of initial geometric imperfection on the nonlinear vibration characteristics of FG-CNTRC rectangular plates excited by a harmonic transverse load. The FG-CNTRC rectangular plate is assumed to be made of the matrix and SWCNTs. The effective material properties of FG-CNTRC plate are estimated by employing the rule of mixture. Hamilton's principle is used to achieve the nonlinear mathematical formulations of governing equations and essential and natural boundary conditions for initially imperfect FG-CNTRC rectangular plates. The provided mathematical formulations are on the basis of the parabolic shear deformation plate theory incorporating the effects of transverse shear deformation and rotary inertia as well as considering the geometric nonlinearity using the von Kármán hypotheses. A multistep solution methodology is utilized to solve the obtained nonlinear partial differential equations in order to provide the frequency-amplitude and forcing-amplitude curves of initially imperfect FG-CNTRC rectangular plates with various edge conditions. In the numerical results, the impacts of initial geometric imperfection, geometries of system, and edge conditions on the nonlinear dynamical responses of imperfect FG-CNTRC rectangular plates are examined. It is concluded that increasing the initial geometric imperfection results in intensifying the softening-type manner of imperfect plate, while no softening treatment can be found in the frequency-amplitude curve of a perfect FG-CNTRC rectangular plate. Also, for the perfect FG-CNTRC plates and imperfect plates with small initial imperfection, the frequency-amplitude curve displays an unstable solution with two limit point bifurcations. However, for imperfect plates with high initial imperfection, it can be found that the nonlinear frequency-amplitude curve displays both softening and hardening behaviors and contains two unstable solutions together with four limit point bifurcations. Furthermore, the number of unstable solution and limit point bifurcation of force-response curve of system strongly depends on the parameters such as the boundary conditions, the plate geometry, and the excitation frequency.

Appendix A

The nonlinear parameters Z_i (i=1, 2, ..., 5) in Eq. (22) are as follows:

(A1)

(A2)

(A3)

(A4)

(A5)

Moreover, the linear elements are as follows:

(A6)

(A7)

(A8)

(A9)

(A10)