Computational analysis of the nonlinear vibrational behavior of perforated plates with initial imperfection using NURBS-based isogeometric approach

In this article, an isogeometric analysis through NURBS basis functions is presented to study the nonlinear vibrational behavior of perforated plates with initial imperfection. In this regard, the governing equations of plate dynamics, as well as the displacement–strain relations, are derived using the Mindlin–Reissner plate theory by considering von Karman nonlinearity. The geometry of the structure is formed by selecting the order of NURBS basis functions and the number of control points according to the physics of the problem. Since similar basis functions are utilized to estimate the accurate geometry and displacement field of the domain, the order of the basic functions and the number of control points are optimized for the proper approximation of the unknown field variables. By utilizing the energy approach and Hamilton principle and discretizing the equations of motion, the vibrational response of the perforated imperfect plate is extracted through an eigenvalue problem. The results of linear vibrations, geometrically nonlinear vibrations, and nonlinear vibrations of imperfect plates are separately validated by considering the previously reported findings, which shows a satisfactory agreement. Thereafter, a coefficient of the first mode shape is considered as the initial imperfection and the vibrational analysis is reexamined. Furthermore, the nonlinear vibrations of the perforated plate with initial imperfection are analysed using an iterative approach. The effects of the perforated hole, initial imperfection, and geometric nonlinearity are also addressed and discussed.


Introduction
Perforated plates are employed in place of their simple counterparts in various industries based on the design criteria and mounting locations of the plate. Thanks to their appropriate strength-to-weight ratio, perforated plates have found extensive applications (Sahu & Dutta, 2002) as filters and purifiers in oil and gas refinery (Baltzer & Norman, 1999;Biliaiev et al., 2015), bridge and ship decks (Kumar & Singh, 2017;Enferadi et al., 2019), sound absorber sheet (Mendez & Eldredge, 2009;Sun & Zhang, 2018), heat exchanger and cooling systems of turbines, engines and nuclear reactors (Dizaji et al., 2018;Wang & Wang, 2020;see Fig. 1). In this regard, these structures are now viewed as one of the most important components in engineering fields. Researchers have applied various numerical, experimental, and analytical methods to investigate the static and vibrational analysis of these plates under different loading and boundary conditions. The first study on the vibrational behavior of perforated plates was accomplished by O'Donnell and Langer (O'Donnell, 1973) where the effective material constants extracted from strain estimation and experimental photo-elastic analysis were used in an "equivalent" solid plate. The work of examining the plate with a square hole was continued with numerical methods (FDM and FEM), analytical approaches (Riley-Ritz), and experimental procedures (Aksu & Ali, 1976;Ali & Atwal, 1980;Reddy, 1982).
Due to the difficulty of analytically solving a plate with a circular hole, the geometry of the plates with square holes was considered equivalently. Afterwards, Burgemeister and Hansen (1996) reported that in the classical equations, the "effective elastic constant" could not exactly determine the resonant frequencies of simply supported perforated plates. By a comparison of two theoretical analyses, a finite element method and a modal test analysis on 13 different perforated plates, they developed a method to calculate the natural frequencies. Grossi et al. (1997) used an optimized Rayleigh-Ritz method to study the vibration analysis of a rectangular plate with a circular hole considering mid-plane rotation. Based on the modal test and FEM analysis, Patil et al. (2007) investigated the free vibration of a rectangular perforated plate with different hole diameters with diagonal pattern. In the modal analysis, fundamental frequencies and their related mode shapes were evaluated by hammer test. The effective natural frequencies were defined by the ratio of natural frequencies of perforated plates with regard to the intact plate. Kwak and Han (2007) presented a new analytical method called Independent Coordinate Coupling to examine the free vibration analysis of a rectangular plate with square and circular holes. In this method, the energies related to the entire plate and hole domain were separately extracted and the best location of holes was determined without changing the sides and hole area ratios. Mali and Singru (2013a, b, c) introduced an analytical model relying on the Rayleigh-Ritz approach for determining the principal natural frequencies of a clamped perforated plate with circular holes. They defined the inhomogeneity in Young's modulus and density caused by equivalent square holes by replacing the circular holes with equivalent square holes using unit step-functions (Mali & Singru, 2013a) and the floor or greatest integer functions (Mali & Singru, 2013b). To approximate the middle surface of the plate, they used a function considering boundary constraints for full clamped conditions. In another research (Mali & Singru, 2013c), they also analysed the linear vibration of rectangular plates with circular holes and rectangular pattern using the energy developed named negative concentrated mass. Holes were determined as the negative concentrated masses in a plate with uniform mass distribution. By applying the boundary conditions and using the cosine approximation function that satisfies the boundary conditions and utilizing the Galerkin method, they succeeded to solve the problem. Rahul et al. (2015) analysed the vibrational characteristics due to different hole shapes of composite plates using the ANSYS software. To determine the stress concentration zone, the finite element analysis was applied and the discrepancy between the simulated results reached 36%. Jhung and Jeong (2015) suggested an equivalent material method for cylindrical and perforated plates and verified their proposed method using the finite element approach on several perforated and intact plates. Considering a nine-node isoparametric element, Kalita and Haldar (2016) investigated the free vibration of a perforated plate under various boundary constraints. Barry and Tanbour (2016) examined the vibration of a perforated plate with rectangular holes with rectangular pattern considering the simple boundary conditions. The natural frequencies were estimated using equivalent properties by means of elasticity theory. In addition, the density and equivalent elastic coefficient were determined by ignoring the changes in Poisson's ratio. Relying on the properties of equivalent material, Luschi and Pieri (2016) followed the previous works and defined a simple analytical model for investigating the vibration behavior of perforated plates being used in mounting sensors and measuring parts. Then, the natural frequencies were calculated for Lame modes and the results were compared with experimental and FEM findings. Hachemi et al. (2017) employed the Mindlin-Reissner plate theory to study the effect of fiber angles, hole diameter, and thickness on the vibration features of a composite plate. Under the influence of fluid flow, Bossi et al. (2017) explored the effect of diameter and number of perforated holes on the fluid-induced instability, mean drag, and lift oscillation coefficients.  presented a theoretical pattern for the sound absorber factor on multilayered elastic micro-perforated plates (MPPs). Lower frequencies and symmetric mode shapes from the vibration of plate as well as sound pressure were calculated under plane wave conditions. The effects of hole diameter, thickness, damping, cavity depth, and perforation ratio were assessed on the sound absorber factor. Moreover, the optimized parameters were suggested to maximize the sound absorber coefficient at a specified frequency band in an MPP with different layers. Konieczny et al. (2019) calculated the stresses in a symmetric, isotropic, and perforated plate with circular holes under the concentrated loads using the Abaqus software and evaluated the results with empirical procedures. Very recently, Sun et al. (2021) evaluated the free vibration characteristics of composite plates with multiple cutouts by employing an FEM analysis with scaled boundary combined with the precise integration method.
Environmental conditions, manufacturing processes, improper handling, or welding can lead to geometrical defects in many structures. The considered geometric imperfections are small deviations of the actual fabricated plate from the intended or designed shape. These imperfections are in the structure geometry and will change the vibration behavior of plates. Hui and Leissa (1983) investigated the initial geometrical imperfection in a large amplitude rectangular plate and showed that an increment in the amplitude of initial imperfection can enhance the natural frequency. Additionally, they proved that if the amplitude of initial imperfection is equal to half of the plate thickness, the hardening behavior will change to a softening one. Using large deflection plate theory, Ostiguy and Sassi (1992) analysed the effect of geometrically initial imperfections on the vibration features of a simply supported rectangular plate exposed to in-plane periodic loads. Han et al. (1994) conducted a nonlinear analysis on multilayer composite plates using the finite element method with a hierarchical approach. After the derivation of the equilibrium equations, the stiffness and tangential stiffness matrices were calculated using the symbolic calculation of high-order polynomial integrals and the Newton-Raphson iteration method. The effect of in-plane displacement on the nonlinear geometrical deformation was also discussed. Ilanko (2002) applied a methodology similar to Hui and Leissa (1983) to deal with the vibration analysis of thin plates with initial imperfection under inplane loading. Using dual Fourier series approximated initial imperfection, he calculated the mode shapes and natural frequencies of the structure. Ribeiro (2003) utilized a hierarchical P-version finite element method to evaluate the nonlinear vibration of thick plates by considering the impacts of geometrical nonlinearity, rotatory inertia, and transverse shear. Moreover, the frequency response was determined utilizing the harmonic balanced method. Amabili (2004) explored the harmonic forced vibration of a rectangular plate at large amplitudes. Considering the initial geometrical imperfection, the experimental results were compared with the numerical data obtained from von Kármán nonlinear plate theory. Subsequently, Amabili (2006) applied the energy method to investigate the nonlinear forced vibrations of a rectangular plate in the presence of mass concentration. The results were illustrated as frequencydomain curves and the time-domain responses. The effect of gravity and initial imperfection was also experimentally addressed. Nonlinear finite element analysis was conducted by Sharma and Kumar (2016) to examine the dynamics of FGM plates with a center hole under in-plane compressive load. In this study, the shear deformation theory and nonlinear kinematic von Karman relation were taken into account. They demonstrated that due to the rigid edge body conditions, the clamped plate with a larger hole could tolerate more buckling load than the plate with smaller holes. Furthermore, it was exhibited that by increasing the hole diameter, the fracture load, maximum shear deformation, and postbuckling stiffness are gradually decreased. Guo and Zhong (2016) introduced a new definition for thin plates and a refined potential energy principle was applied to deploy the improved imperfection in static experimental measurements. Ghayesh (2019) inspected the nonlinear behavior of an imperfect microscale Timoshenko beam made of functionally graded (FG) material using the Kelvin-Voigt viscoelastic model for internal viscosity. By extracting the stress-strain equations using the Hamiltonian principle, three sets of coupled equations of motion were achieved. In their study, the Galerkin procedure was employed and the resultant nonlinear equations were analysed by applying the pseudo-arclength iterative method. By interpreting the frequency diagrams, it was observed that large geometrical imperfection will lead to initial softening followed by a hardening behavior. Ruocco and Mallardo (2019) combined the finite element and finite strip methods to evaluate the buckling and vibration characteristics of nanoplate with initial imperfection relying on the Mindlin-Reissner plate theory. For a perfect description of the geometrical imperfection, they applied the finite element method for incorporating the out-of-plane and in-plane buckling equations. They also deployed the Kantorovich method to solve nonlinear equations. Several numerical examples were addressed to illustrate the sensitivity of the nonlocal parameters with regard to geometrical imperfection. Hassani and Javanbakht (2021) proposed the modified control method to assess the failure process of a perforated panel with a web-reduced beam section considering various geometrical modifications through finite element analysis. Foremost, it should be also mentioned that the closed-form solutions for free vibration analysis of Mindlin plate are complicated and restricted to special cases because of the number of governing equations, independent coordinates, and kinetic parameters involved (Xing & Liu, 2009).
To this end, various numerical methods are used by investigators to study free vibration analysis of plates and beams including residual power series method (Jena et al., 2021), Navier's approach (Jena et al., 2020), differential quadrature element method (Liu & Liew, 1998;Liu et al., 2017), finite element method and boundary element method (Nguyen-Xuan et al., 2008;Chau-Dinh & Le-Tran, 2020;Fallahi et al., 2020;Torky & Rashed, 2020;Vo & Ton-That, 2020;Yang et al., 2020), and radial point interpolation method (Cui et al., 2010). Nowadays, meshfree methods are developed, enriched, and applied to analyse a wide variety of problems with complicated geometries (Jalal et al., 2019;Hosseini et al., 2021). In addition, Guo et al. (2019) examined the bending analysis of a thin plate with various geometries using a new method known as the deep collocation method through a deep neural network (DNN) and computational graphs. However, there are still many challenges for the DNN-based methods including the effect of different types  (Ostiguy & Sassi, 1992). of neural networks, activation functions, and in the optimization problems, the analysis of approximation spaces is so difficult that even linear problems become nonlinear and nonconvex discrete problems Samaniego et al., 2020).
An in-depth review on the perforated plate and Mindlin plate literature reveals that the analytical and numerical analysis of perforated Mindlin plates with different shaped holes is a complicated and time-consuming procedure, which will be even harder by considering the initial imperfection, as a geometrical nonlinear case. In order to deal with the complex structures and geometries, along with the mesh-free and DNN-based methods, the isogeometric analysis was presented by Hughes et al. (2005) in 2005. NURBS functions are very efficient to implement and well used in computer-aided geometric design. Furthermore, these functions can naturally deal with the uncertainties and optimization problems within a single framework (Schillinger et al., 2012;Bu¨nger et al., 2020). Regarding this, recently, a lot of works have utilized NURBS-based IGA to investigate the static bending, vibration, and buckling analysis of different geometries and structures based on various theories of shear deformation, loading, and boundary conditions that can be similar or extended to recent literature (Thai et al., , 2015Allam et al., 2020;Bekkaye et al., 2020;Arshid et al., 2021;Bakoura et al., 2021;Bellifa et al., 2021;Guellil et al., 2021). Nguyen-Thanh et al., (2015) investigated a cracked thin shell based on the Kirchhoff-Love theory using extended isogeometric analysis. Chen et al. (2017) analysed the in-plane vibrations of various orthotropic plates including a rectangular plate with circular and square holes, rhombus, and trapezoidal plates.  studied the natural frequencies of bidirectional FG beams; Liu et al. (2019) investigated the static bending and buckling analysis of FG plates using a refined plate theory (RPT); Fan et al., (2020) examined the bending analysis of laminated composite and carbon nanotube-reinforced composite based on higher order shear deformation theory; Singh & Bhar (2020) investigated the porosity-dependent nonlinear behavior of a porous FG rectangular micro-plate; and Zhong et al. (2021) analysed the free vibration and bending of FGM plates with various cutouts considering the thermal effects by defining IGA on high-order splines. Cuong-Le et al. (2021) proposed a 3D numerical weak form approach to deal with the vibration and buckling analysis of FG plates and shells based on the four-variable RPT; Phung- Van et al. (2021) analysed the vibration characteristics of multilayer FG nanoplates. Although NURBS basis functions have some limitations for representing the complex and multiply connected domains and performing on local refinement due to its tensor-product forms, however, multiple patches give more flexibility in the construction of domains and make local refinement more practicable depending on how patch interfaces are jointed. In a multipatch geometry, distinct patches are modeled the complex body and keeping the shape of the elements inside the patches (Hughes et al., 2005). In addition, a new approach for local order elevation as well as weak enforcement of continuity was applied to multipatch with nonconforming interface (Fazilati & Khalafi, 2021;Yuan et al., 2021). To attach different NURBS patches, Zhao et al. (2017) analysed the plate made up of multiple NURBS patches with nonconforming interfaces by employing the Nitsche method.
According to the aforementioned literature survey and to the best authors' knowledge, there is no essential effort to analyse the nonlinear vibrational behavior of perforated Mindlin plates with initial imperfection in the context of IGA, and therefore, this study is aimed at exploring the superiority of this efficient numerical approach to deal with the simultaneous effects of complicated geometries (i.e. perforation) along with the realistic nonlinearity (i.e. initial imperfection). To demonstrate the effectiveness of the present method, several numerical studies are conducted and the influence of hole size, perforation pattern, geometric nonlinearity, and initial imperfection is thoroughly addressed and discussed.

Theory
The studied perforated plate is considered to be isotropic with the respective length, width, and uniform thickness of a, b, and h (Fig. 2). The transverse shear deformation and rotary inertia based on Reissner-Mindlin's formulation are also included in our assumption (Fung, 1977;Chakraverty, 2008;Draiche et al., 2019).
It is supposed that the displacement vectors of u and v at each point of the plate domain in the x and y directions can be expressed in terms of the mid-plane translations based on FSDT. However, depending on the physical domain and boundary conditions of the problems, different shear deformation theories are used Bousahla et al., 2020;Matouk et al., 2020). Hence, φ y and φ x are remarked as the independent rotations of the normal to the neutral surface around the y and x axes (see also Fig. 2). The displacement vector along the z-direction, shown by w, is related to the initial geometrical imperfection w I . These initial imperfections of plate elements are normal displacement dependent on zero initial stress and therefore the in-plane initial imperfections are neglected. Thereby, the displacement field of FSDT is obtained by setting c = 0, in equation (1), considering the geometric function together with the distribution of transverse shear strains and stresses through the plate thickness (Zienkiewicz et al., 2014;Draiche et al., 2019;Bourada et al., 2020;Bousahla et al., 2020;Matouk et al., 2020;Rouabhia et al., 2020): where u 1 , u 2 , and u 3 are described as the functions of (x, y, z, t), the displacement of mid-plane, and φ x , φ y are not dependent on z.
The strain-displacement relations can be described by By substituting equation (1) into the strain-displacement equations, the strain relations can be expressed as where the plate curvatures can be defined by Relying on the von Karman assumption, the derivations of u and v with regard to x, y, and z are small; hence, w is independent of the z-direction and Green's strain can be rewritten based on mid-plane deformation. Therefore, the strain-displacement equations can be rewritten as (Amabili, 2008;Zienkiewicz et al., 2014;Rouabhia et al., 2020) where and The stress constitutive equations for the isotropic element of the plate are then given by (Zienkiewicz et al., 2014) For the plane stress case, one obtains In equations (7) and (8), the plate element coefficients in [Ĉ] are the diminished stiffnesses that depend on the Poisson's ratio, ν, elastic modulus E, and the shear modulus, G. Moreover, the subscripts 1 and 2 represent the principal directions of the isotropic plate in the x and y directions, respectively. One can indicate that (Piegl & Tiller, 1996) ⎧ where the bending strain is described as Thereby, the in-plane stresses and moments of M can be calculated by equation (13): and the in-plane shear stresses Q x and Q y are extracted as in which the membrane, flexural, and coupling rigidities can be determined as (Zienkiewicz et al., 2014) In this research, only the symmetric isotropic plate is considered, and therefore, according to equation (15), the elements of the matrix [B] are assumed to be zero and in-plane coupling between the transverse bending and stretching is no longer valid. Moreover, the matrices [D] and [A] are given by (Chakraverty, 2008) [ in which i, j = 1, 2, 6. To establish the governing equation in weak form, Hamilton's principle is implemented as ) where δ represents the variational operator. I , S , and K are potential energy from initial stress, strain, and kinetic energies, respectively. For each plate element, the total potential energy can be written as The variation of strain energy δ S is also calculated as while the variation of kinetic energy δ K is determined as Due to the absence of any initial stress, one has Substituting equations (20)-(22) into equation (6) results in the following dynamic equilibrium equation in the weak form: Using the generalized coordinates, we can also reduce the coordinates of u and v directions, which will be explained in the next section (Zienkiewicz et al., 2014): in whichm Matrix summation is performed with the corresponding displacement or rotation. Up to this point, the equations are extracted based on the intact rectangular plate element. To demonstrate the effect of holes in the perforated plate on the obtained relations, at first, the isogeometric method is defined and then the necessary corrections are made by mapping the rectangular element to the desired element in the perforated plate.

Isogeometric analysis approach
Similar to the isoparametric approach, the same NURBS base functions are used in the isogeometric method to analyse the model and approximate the physical geometry of the body and the field variables . The parametric space can be defined by one or more patches, and each patch itself is considered as a subdomain or macro element. B-spline basic functions are defined in a one-dimensional parametric domain using a set of parametric coordinates of nondecreasing real numbers called knot vector in the form of H = {η 1 , η 2 , η 3 , . . . , η n+k+1 }, in which k is the polynomial of B-spline, η j shows the jth knot, and n represents the number of control points or base functions. The interval between two consecutive knots is called either the knot span or element and as the knot span is similar, the knot vector is uniform and in other cases, it would be a nonuniform vector. Each knot vector contains 'n−k' knot span that will be zero as the values are repeated. Having a knot vector and the order of basic functions, one-dimensional B-spline basic functions are determined by the Cox-de Boor recursion formula as follows (Piegl & Tiller, 1996;Cottrell et al., 2009): For k = 1, 2, 3, ..., we have If the denominator of the fraction of some coefficients in equation (27) is equal to 0, the values of these coefficients are also considered as 0. The jth derivative of a k-order polynomial can be represented as Differentiating both sides of equation (28) with respect to ξ results in the higher order derivatives as follows: Expanding equation (29) by equation (28) leads to the extraction of lower order relationships. The knot coordinates are the same as the coefficients of the basic functions. Using the basis function of order "k" and the associated control points or position vectors, it is possible to create a piecewise-polynomial B-spline curve in the physical coordinate space via the following linear combination: The control points corresponding to the basic functions characterize the geometry. The control polygon is then formed by connecting the control points with piecewise linear interpolation lines. NURBS basis functions are derived from the generalization of B-splines in which B-spline is given in R d+1 , mapped onto a reduced space R d , in which "d" represents the space dimension. NURBS basis function R k j (η) can be subtitled in a weighted basic function: in which w l is the l-th weight coefficient, 0 < w l ≤ 1, and W(η) shows the weighting function. Having projective curve C w (η) and its associated projective control points CP w l , the control points for the NURBS curve are derived by wherein (C P l ) k is the kth membership of C P l . Rational surfaces can be constructed in terms of the rational basic functions by Applying the tensor product of two NURBS basis functions by p and q orders in ξ and η directions as N k, r (ξ ) and M l, s (η), k = 1, . . . , r and l = 1, . . . , s, the NURBS surface via r × s net of control points CP k, l is computed as NURBS basis functions have some limitations for representing complex, multiply connected domains and performing local refinement due to its tensor-product forms (Chen et al., 2017). However, multiple NURBS patches are necessary to describe the domains. The Analysis may frequently demand higher orders than geometric design; however, in general, the tendency is to utilize the least polynomial order in each parametric direction. The control points corresponding to the basic functions construct the geometry. To replicate the circular features of each plate element with a hole, set r = 2 and s = 2, although there is no ambiguity in our selection (see Fig. 3), the knot vector selection is carried out by determining the necessary elements and the required level of continuity across each element and patch boundaries.
For the first attempt, the knot vectors are used as H = = {0, 0, 0, 0.5, 1, 1, 1 } . It is worth noting that the "pre-image" of the NURBS projection is produced by considering only the nonzero knot spans (Fig. 3a). This physical space is a patch with 16 control points (Fig. 3b). In order to create a square plate element by four patches (Fig. 4), three spaces are defined in the isogeometric method namely the index, parametric, and physical spaces. Figure 3b indicates these three spaces along with their corresponding mapping. The resulting element is known as the knot span in the parametric domain, i.e. [ξ k , ξ k+1 ] × [η l , η l+1 ]. Figure 3b illustrates the mappings of S :ˆ → together withφ :˜ →ˆ from parametric to physical space and from parent to parametric space. Thus, the physical space can be mapped onto the parent element using the composition of the inverses of two mappings, i.e. X = S •φ. The mapping ofˆ e = [ξ k , ξ k+1 ] ⊗ [η l , η l+1 ] for nonzero parametric space is also presented from parent to parametric space (Samaniego et al., 2020) as follows:φ The corresponding Jacobian matrix determinant can be subtitled in The transformation Jacobian for the mapping of S :ˆ → yields Then, the Jacobian can be determined by in which N is the total number of the basic functions. Accordingly, the last Jacobian for the combined projection of X e = S •φ is given by By the combinatorial projection and utilizing the last Jacobian, the conversion equation for the expression of Gaussian rule between the physical spaces and parent element can be derived as in which N el stands for the quantity of nonzero elements (Cottrell et al., 2009). The integral outcome is derived using the conversion formula by replacing the variables and coordinates of control points into the virtual work expressions. There have been many investigations conducted on extracting the optimal quadrature method to numerically construct the mass and stiffness matrices, even though there are no general quadrature rules for spline spaces of arbitrary order and continuity to integrate them Bartoň et al., 2017;Hiemstra et al., 2017). Hence, an algorithm is exploited that computes the optimal and refined quadrature rules for NURBS basis integrands that can be applied for any uniform and nonuniform knot vector of necessary order and continuity (Barton & Calo, 2016;Johannessen, 2017).

Refinements
In the isogeometric analysis, once the initial physical mesh is constructed, the actual geometry in physical space is formed. The hrefinement, p-refinement, or k-refinement procedures can be employed to improve the mapping methods. At each refinement step, there is no need to form a geometry again. In the h-refinement method, similar to the FEM method, a new curve or surface is created by inserting a new knot into a knot vector, which is parametrically and geometrically similar to the primary curve; however, its control point coordinates are moved. The p-refinement method is performed by enhancing the polynomial order of the basic functions to demonstrate the geometry. To maintain the C 0 continuity in the boundaries and C p−k in the internal knots, without inserting new knot values, it is necessary to insert each knot once again within the knot vector. In addition, the shape geometry in the parametric space does not change in the process of knot insertion. The method of order elevation starts with repeating all knots within the knot vector, which continues until their iteration equals the order of the polynomial. The curve is then divided by adding extra knots, increasing their polynomial order, and finally removing the extra points to form the enhanced order in the B-spline curve. If a new knot value ξ i is placed between two consecutive knot values of B-splines curve of order p, the continuity of the pth derivatives of the basic functions in the new knot will be p−1. In order to satisfy a discontinuity in that knot while having order raising in the (p−1)th derivative, the value of each knot should be repeated as long as the continuity of the derivatives remains (p−1). If the value of the knot ξ i is added, the base functions at point ξ i will possess a continuum of derivatives of order q−1; this process is called k-refinement (Piegl & Tiller, 1996). NURBS has also shown significant improvement in the analysis of structural vibration problems in terms of robustness and accuracy, by using k-refinement as compared to higher order p-refinement . IGA was initiated to rely on NURBS (Piegl & Tiller, 1996), and then, in order to develop the local refinement, has been outstretched to T-splines (Bazilevs et al., 2010), hierarchical approach , and polynomial/rational splines over hierarchical T-meshes (PHT/RHT-splines; , 2017Lu et al., 2018). After that, Anitescu et al. (2018) defined an adaptive local refinement using a recovery-based error estimator that relies on PHT/RHT-splines to resolve the singularities that may appear while the refinement propagates throughout the computational domain. Although the continuity is less than the highest amount of CP-1 generated by NURBS, the diminished smoothness facilitates the fabrication and refinement processes and efficiently reduces the computational time for direct and iterative solvers (Allen & Al-Qarra, 1987). Multiple patches give more flexibility in the construction of the domain and make the local refinement more practicable depending on how patch interfaces are joined. If the control points at the interface are in the same location in the space, both sides of the interface domain between two patches require the same NURBS basis functions to satisfy the continuity condition and mesh refinement from one patch to another (Fig. 4b). Another possibility may be the situation in which the number of control points of each patch at the common interface is different, whereas it would be a knot-refined version of another one (Fig. 4c). In this situation, C0-continuity is satisfied in the solution across the interface domain and higher continuity will also be carried out by utilizing similar involvement equations in the normal direction.

Weak form of vibration equations
The mid-plane displacements and rotations for each element of the plate in generalized coordinate can be described by 2D-NURBS basis functions as (Zienkiewicz et al., 2014) In equation (41), the generalized coordinate vector p is defined by { p} = {p u , p v , p w, p φy , p φx } T and R is a diagonal matrix of NURBS basis functions as described in equation (34); each of these elements is proportional to one direction of the displacement vector u with respect to the generalized coordinates. So, in 2D space, R(x(ξ, η), y(ξ, η)) is defined as in which {R u } = {R u1 , .., R uj , . . . , R unCP } T , wherein nCP stands for the number of control points in each knot span; then, R u = R v are similar to the rest. The element's stiffness, mass, and force matrices are constructed by replacing the displacement field estimation into the weak form. The linear in-plane bending, initial imperfection, and transverse shear strains are displayed in equation (43) as where Therefore, the in-plane membrane, bending, and shear stiffness matrices are illustrated through equations (45) to (47): and the mass matrix is given by where [C M ] is defined in equation (21). Substituting equations (41)-(48) into equation (23) leads to In the linear stiffness matrix, the expressions of K L 13 and K L 23 appear due to the initial imperfection, and the expression for K L 33 is normally available in the plate equations. They are now influenced by the initial imperfection. In the nonlinear stiffness matrix, the expression for K NL33 is linearly affected by the initial imperfection. The expressions of mass and stiffness matrices are presented in Appendix A. For the cases where the transverse displacement is much higher than the mid-plane and in-plane displacements, it is reasonable to ignore the mid-plane inertia, for the sake of simplicity. Consequently, the generalized in-plane displacements can be simplified as follows: The expressions of K L S33 and K NL S33 are affected by nonlinear expressions and initial imperfection expressions as explained before; the effect of each of the mentioned expressions on the stiffness matrix can be calculated (see Appendix A). Finally, the matrix form of the equation of motion is written by (Lau et al., 1984;Allen & Al-Qarra, 1987;Amabili, 2008): To numerically solve the abovementioned equations, the nonlinear stiffness and initial imperfection terms must first be set to zero, and then linear free vibration analysis be carried out (equation 52b).
−ω 2 Subsequently, the initial imperfection and nonlinear coefficients are estimated to deal with the nonlinear governing equations using an iterative method. However, the forced vibration equation subject to F = q 0 cosωt can be rewritten as follows (Allen & Al-Qarra, 1987;Amabili, 2008): in which q =[q 0 0 0] T cosωt.

Newton-Raphson numerical scheme
The nonlinear eigenvalue problem can be solved through an iterative method (equation 52). In this regard, after forming the stiffness and mass matrices, the problem is solved via computational software by considering relevant boundary conditions, while ignoring the initial imperfection and nonlinear terms. Afterwards, a coefficient of the first mode shape is considered as the initial imperfection and the problem is solved again by neglecting the nonlinear terms. After this step, the nonlinear vibrations are analysed by applying the Newton-Raphson method via an iterative procedure. For a nonlinear system, the element equilibrium equation would be (Lau et al., 1984;Collier et al., 2013) in which "Res" is residual and F denotes the generalized forces stem from the variation of applied external forces. The solution process for the assembled nonlinear equation (52) is based on the Newton-Raphson method that is composed of a series of linear equations. Thus, the Taylor series expansion of residual Res(p i+1 ) in the neighbor of p i is while the subscript T represents the tangent in K T stiffness as Substituting equation (53) into equation (55) results in . Then, the tangent stiffness matrix K T is Considering the initial imperfection in both linear and nonlinear stiffness terms of equation (51) referring to Appendix A, the relation for K T is given by Subsequently, the boundary conditions are applied and equations are simulated by an iterative Newton-Raphson method. It is worth mentioning that in this procedure, the nonlinear terms are linearized around the equilibrium points by The first-order approximation of residual function yields The series of consecutive approximation is presented by Here, n + 1 is the number of equilibrium solutions and i represents the Newton iteration sequence. The iterations at each process are followed until the residual becomes less than the desired tolerance.

Results and Discussion
In order to validate the results of the present analysis, the linear dynamic behavior of a partially perforated plate is examined and the first nine natural frequencies are tabulated in Table 1 for different hole radii and then the fundamental frequencies are compared with the consequences of Mali and Singru (2013c) obtained from numerical and analytical methods, as listed in Table 2.     Table 2, it is revealed that the maximum error value is 1.04%, which is satisfactory reflecting the soundness of the current analysis in determining the vibrational features of the perforated plate. As illustrated in Fig. 5, the relative error of the present analysis is much less than the other methods and is in excellent agreement with FEM results in comparison with the other methods. It is also observed that with increasing the diameter of the perforation hole the fundamental frequency shifts upwards accordingly. To further justify the results of the current approach, the linear dynamic behavior of a plate with initial imperfection is simulated and compared. For this purpose, the results of this work are evaluated with Amabili (2006) and tabulated in Table 3. It is confirmed that the greatest relative error is below 1%, demonstrating the precise modeling and analysis in this research. In this study, the nonlinear governing equations are solved by an iterative algorithm on the basis of the Newton-Raphson method. In the final step of verification of the results, the nonlinear vibrational analysis of an isotropic plate is conducted using the proposed method after four iterations and the outcomes are compared with those reported in Lau et al. (1984), as given in Table 4. It is exhibited that the greatest relative error is 0.79%, confirming the excellent accuracy of our analysis in determining the nonlinear behavior of the considered structure. Subsequently, the nonlinear vibrations of the perforated plate with initial imperfection are investigated. The specifications of the perforated plate considered in the simulations are listed in Table 5. The schematic configuration of the perforated plate with circular holes with triangular and square patterns is depicted in Fig. 6. To show the efficiency of the isogeometric analysis method in analysing the vibrations of perforated plates, at first, a plate with a central hole is modeled using four patches including a quadratic knot vector due to the second degree of the circle shape with 24 control points (Fig. 4).   The shape geometry is precisely modeled by defining the initial knot vector with the minimum number of control points and the order of corresponding basic functions (see Fig. 4a). To determine the appropriate number of optimal control points and the order of basic functions to attain acceptable accuracy at a smaller computational burden, the necessary studies are conducted. In this research, by preparing a Matlab code, the number of control points and the order of basic functions are increased separately and simultaneously (Fig. 7). In the current analysis, the number of optimal control points together with the appropriate order of the basic functions can be determined to justify the convergence and accuracy of our findings for vibration analysis of perforated plates. These steps are shown in Fig. 7 by performing the convergence of h, p, and k refinement. Upon uniform distribution of the holes (instead of concentration at the center), the influence of hole size on the vibrational characteristics is explored (as indicated in Fig. 6) and the results of the first 12 natural frequencies are listed in Table 6. According to the presented results in Table 6, it is found that an increase in the diameter of uniformly perforated holes results in reducing the natural frequency. As demonstrated in Table 2 and Fig. 5, the frequency tends to significantly increase due to the appearance of the perforation holes in the center of the plate as the diameter becomes larger. It should be noted that for modes 5, 6, 9, and 10, by increasing the diameter of the perforation holes to more than 16 cm and consequently increasing the perforated surface, it behaves like a plate with a central hole again and the frequency is increased. This is due to the fact that the stiffness of the plate is affected by the location and the number of holes on the side in this form of higher order modes (see Fig. 8).
The effect of the triangular and rectangular patterns under clamped boundary conditions is presented in Table 7. The mode shapes for the clamped perforated plate are illustrated in Fig. 8. In addition, according to the presented results in Table 7, it is inferred that the perforation patterns with uniform distribution of perforation holes do not have much effect on the first natural frequencies, but for simply supported conditions, the effect of rectangular pattern on the reduction of 10th frequency is up to 8.12% more than the triangular pattern. Afterwards, the effect of initial imperfection on the variation of the natural frequencies of perforated plate is also examined. A new parameter containing the amplitude of the first mode with respect to the thickness of the plate is defined to analyse the influence of the initial imperfection. The simulated results for the natural frequencies of perforated plates with different initial imperfections are presented in Tables 8 and 9 for clamped and simply supported boundary conditions. One concludes that when the amplitude ratio takes the larger values, the nonlinear frequency due to initial imperfection increases monotonically. By  increasing the similar amplitude ratio for simply supported and clamed perforated plates, the frequencies percentage variation is up to 6 times higher. This is due to less maximum displacement of clamped perforated plate compared to simply supported one. The variation of the nonlinear to linear natural frequencies of both simply supported and clamped perforated plate with respect to initial imperfection is displayed in Fig. 9 comparing the results of flat plate and perforated one. The relative frequency variations of both simply supported and clamped perforated plate are greater than those of the intact plate, although these variations are less in the clamped one. Moreover, Fig. 10 indicates the effect of thickness ratio (a/h) on the variation of the first 10 natural frequencies of the perforated plate. According to the illustrated results in this figure, it is concluded that as the thickness ratio takes the higher values, the natural frequencies exponentially decay. In this example, the effect of initial imperfection on the ratio of nonlinear to linear frequencies is investigated by varying the thickness ratio, which is not shown due to its negligible effects. In the next case study, the effect of initial amplitude (w max /h) on the nonlinear natural frequencies of the system is studied. To this end, the ratios of the nonlinear to linear natural frequencies are listed in Tables 10 and 11 for clamped and simply supported perforated plates, respectively. As presented in Tables 10 and 11, the natural frequencies are enhanced by increasing the initial amplitude of vibration, which is higher for a simply supported case compared to clamped boundary conditions. For more clarification, the variations of the nonlinear frequency ratios of the clamped and simply supported perforated plate are graphically depicted in Fig. 11 as a function of the initial amplitude ratio. As described in Tables 9 and 10, the ratio of higher modes to the first frequency and the nonlinear variation of frequencies in the simple support is higher than the clamp. The hole perforation also reduces the stiffness   and increases the effect of geometric nonlinearity. Moreover, the convergence diagram of the iterative method for the nonlinear solution is provided in Fig. 12. It is demonstrated that after eight iterations, a satisfactory convergence is achieved and the relative error approaches to 0 after a small number of iterations confirming the efficiency of the proposed approach. To further explore the nonlinear vibrational characteristic of a perforated plate with initial imperfection, the nonlinear to linear natural frequencies of perforated plates with initial imperfection are also investigated. For this purpose, Table 12 shows the first 10 nonlinear frequencies of the imperfect perforated plate with simply supported conditions in which the initial imperfection is assumed to be similar to the first mode       shape and takes a maximum amplitude of 0.2 with respect to plate thickness. The results of a clamped perforated plate are also presented in Table 13. One can infer that when the initial amplitude becomes larger the nonlinear frequency ratio is also increased and the variations are more remarkable in the case of simply supported boundary conditions. Afterwards, the effect of forced vibration with harmonic forcing function on the variation of vibration amplitude of perforated plate is performed using frequency response diagrams for different perforation hole diameters and initial imperfections, as shown in Fig. 13a and b. As indicated in Fig. 13, for the perforated plate with simply supported and clamped boundary conditions, by increasing the amplitude of the harmonic forcing function, the amplitude of vibration increases; however, the hardening behavior of the structure decreases, although the former is greater than the latter. In addition, the effect of perforation hole diameter is to decrease the hardening features of perforated plates (Fig. 13b). Finally, the first nonlinear frequency of the simply supported perforated plate with different initial imperfections with respect to vibration amplitude ratio is plotted in Fig. 14a and b for the free and forced nonlinear vibrations, respectively. As depicted for   the free vibrations in Fig. 14a, it is revealed that initial imperfection makes the positive (outward) and negative (inward) amplitudes asymmetric, in contrast to the intact plate. Also, the nonlinear natural frequencies are reduced showing the so-called softening type, although the nonlinearity seems to be a hardening type (i.e. the amount of vibration increases as the amplitude increases) when the amplitude ratio of vibration changes in a certain interval between 0.4 and 1. Furthermore, as shown in Fig. 14b, when the values of initial imperfection increase, some different phenomena are found in the dynamics of the system. When there is no initial imperfection (w I = 0), or at lower values of initial imperfection, the system shows robust hardening behavior in which the hardening effect is remarkable by decreasing the amount of imperfection. On the other hand, by further increasing the initial imperfection, the system may exhibit softening behavior at some intervals of initial amplitude and become stiffer at some other intervals of initial amplitude.

Conclusions
In this paper, the nonlinear vibrations of the Reissner-Mindlin perforated plate with initial imperfection were analysed considering von Karman nonlinear approach and the NURBS basis function isogeometric method. The order of the NURBS basis functions and the optimal number of control points were determined to obtain the desired geometry and approximate the displacement field with high efficiency and fast convergence. The geometric nonlinearity and initial imperfection terms were applied to the displacement relations and the equations of motion. The eigenvalue equations with nonlinear stiffness were established using the energy method and Hamilton principle. After full validation of the results with the existing models, several studies related to perforated plates, imperfection plates, and geometrical nonlinear plates were investigated using the present method. The influence of perforation hole size and pattern, geometric nonlinearity, and initial imperfections on the nonlinear vibrations of the perforated plate was explored.
The results indicated that an increment in the diameter of the perforated plate reduces the natural frequencies in the case of the perforated plate with distributed holes, although a relative increase in the perforated area section of the plate to more than one-third of the total due to the shape of the mode, for some high-order mode shapes, causes the plate to behave similarly to a plate with a hole and increased the frequencies in the plates. Any decrease in the thickness causes the decrement in the frequencies too but its effect on initial imperfection variations is negligible. Any increase in the vibration amplitude enhanced the influence of the geometric nonlinearity and initial imperfection on the natural frequencies of a perforated plate and this effect was more considerable for simply supported plates compared to clamped supports. By increasing the amplitude ratio, the frequency percentage variations of simply supported perforated plates were up to 6 times higher than clamped ones. When the initial amplitude became larger, the nonlinear frequency ratio was also increased and the variations were more remarkable in the case of simply supported boundary conditions. The hardening behavior of the perforated plate changes by the increasing amplitude of harmonic forcing function, perforation hole diameter, and initial imperfection. Free nonlinear vibrations of the perforated plate additionally altered with increasing the initial curvature. Despite the increase in the linear frequencies, the symmetry and equality of nonlinear frequencies vanished for similar outward and inward amplitude ratios. When the amplitude of vibration increases, its initial softening-type behavior may change to the hardening type. 0 0 0 The linear stiffness matrix is decomposed based on the effects of in-plane forces and initial imperfection and equations (45)-(47) as follows: Expressions that are affected by the initial imperfection include K 13 L 3 , K 23 L 3 , and K 33 L 3 , which are where The nonlinear expression is The nonlinear expression K 33 NL2 , in addition to being affected by the initial imperfection, is also linearly dependent on the transverse curvature: Downloaded from https://academic.oup.com/jcde/article/8/5/1307/6370800 by guest on 21 September 2021 The expression K 33 L S is K 33 L S = K 33 L 1 + K 33 L 2 + K 33 in which is also linearly and squarely dependent on transverse displacement and can be defined as where