Static and dynamic analysis of functionally graded fluid-infiltrated porous skew and elliptical nanoplates using an isogeometric approach

Abstract

The analysis of static bending and free and forced vibration responses of functionally graded fluid-infiltrated porous (FGFP) skew and elliptical nanoplates placed on Pasternak’s two-parameter elastic foundation is performed for the first time using isogeometric analysis (IGA) based on the non-uniform rational B-splines (NURBSs) basis function. Three types of porosity distributions affect the mechanical characteristics of materials: symmetric distribution, upper asymmetric distribution, and lower asymmetric distribution. The stress–strain relationship for Biot porous materials was determined using the elastic theory. The general equations of motion of the nanoplates were established using the four-unknown shear deformation plate theory in conjunction with the nonlocal elastic theory and Hamilton’s principle. A computer program that uses IGA to determine the static bending and free and forced vibration of a nanoplate was developed on MATLAB software platform. The accuracy of the computational program was validated via numerical comparison with confidence assertions. This set of programs presents the influence of the following parameters on the static bending and free and forced vibrations of nanoplates: porosity distribution law, porosity coefficient and geometrical parameters, elastic foundation, deviation angle, nonlocal coefficient, different boundary conditions, and Skempton coefficients. The numerical findings demonstrated the uniqueness of the FGFP plate’s behavior when the porosities are saturated with liquid compared with the case without liquid. The findings of this study have significant implications for engineers involved in the design and fabrication of the aforementioned type of structures. Furthermore, this can form the basis for future research on the mechanical responses of the structures.

Graphical abstract

Keywords

static bending / free and forced vibrations / nonlocal theory / isogeometric analysis / fluid-infiltrated porous nanoplates

1 Introduction

Isogeometric analysis (IGA) is a computational method that provides the opportunity to include finite element analysis (FEA) in traditional non-uniform rational B-splines (NURBSs)-based CAD design tools. This is possible via IGA. To assess new designs during the development phase, the data must be translated between CAD and FEA software. This is a difficult operation because the two computational geometric techniques are distinct. IGA directly incorporates a sophisticated NURBS geometry into FEA application. This is the foundation of the majority of CAD products. Hence, the use of a single dataset is possible throughout the process of designing, validating, and refining of models [1–3]. Within the FEM framework, a significant number of studies extensively used B-splines as trial and test functions. IGA is the term used for the unification and presentation of the framework, and it was developed by Hughes et al. [4,5]. Isogeometric techniques are widely used in computational methods, particularly for modeling the mechanical behavior of complex and intelligent systems. Cottrell et al. [6] shared personal Argyris memory. The sequel revives the geometrical spirit of Argyris by applying IGA to difficulties in structural vibrations. After introducing the IGA, it was used for rods, thin beams, membranes, and thin plates. This approach uses three-dimensional solid models, rotationless beams, and plate models. Hughes et al. [7] addressed NURBS-based isogeometric quadrature rules. The “half-point rule” asserts that optimum rules include a number of points equal to half the number of degrees of freedom or basis functions of the space being investigated. The “half-point rule” describes this guideline. Polynomial order does not affect the half-point rule. The smoothness of the basis functions across the element boundaries must be considered for effective rules. They used numbers to determine efficient and practical rules. Dörfel et al. [8] illustrated the contribution of T-splines. Specifically, T-splines allow T-junctions, which are dangling nodes in the FEM. Rectangular patches in the T-spline parameter space can be split using simple rules, allowing local refinement while maintaining the geometry. Modern posteriori error estimation can be coupled with T-spline refinement. Buffa et al. [9] discretized the two-dimensional Maxwell equations. Inspired by the novel paradigm of IGA, Hughes et al. [4] suggested an approach based on bivariate B-spline spaces that are well-suited for electromagnetics. They generated parametric B-spline spaces with varying inter-element regularities. They used these spaces (together with physical push-forwards) to solve Maxwell source and the eigenproblem. In 2010, Bazilevs and Akkerman [10] employed the residual-based variational multiscale (RBVMS) to calculate Taylor–Couette flow at a high Reynolds number. They showed that the RBVMS formulation conserves angular momentum globally, which is critical for rotating flows but not shared by traditional stabilized formulations. Weak Dirichlet boundary conditions increase the RBVMS accuracy near solid walls with shallow turbulent boundary layers. Conservative boundary forces and torques are calculated for border situations with poor enforcement. Spatial discretization uses NURBS-based IGA, and mesh refinement evaluates convergence. Cohen et al. [11] explained that an equivalent notion of model quality exists in IGA and that this concept behaves similarly to mesh quality, which is used in traditional FEA to characterize the effect of the mesh on the analysis. As a consequence of these findings, “analysis-aware modeling” must be added to the discipline of modeling to ensure that IGA is more accessible. Valizadeh et al. [12] examined the static and dynamic properties of functionally graded material (FGM) plates using a non-uniform rational B-spline-based iso-geometric finite element approach. Dsouza et al. [13] introduced a non-intrusive technique paired with a nonuniform rational B-spline-based isogeometric finite element method to investigate the stochastic static bending and free vibration characteristics of FGM plates with innate material unpredictability. Hu et al. [14] provided a local B-bar formulation that addressed locking in a degraded Reissner–Mindlin shell formulation within the framework of IGA. Specifically, IGA has been utilized in several areas [15–23].

In the realm of computational mechanics, mechanics has gained considerable attention in the application of IGA methods to assess the mechanical response of beams, plates, and shell structures, which are constructed of sophisticated materials. Since the suggestion of using this approach, there has been a significant increase in the number of scientists engaged in research and development, especially in the middle of the third decade of this century. Specifically, micro- and nano-structures have been extensively researched and used in many different aspects of daily life. The following are noteworthy achievements. Ansari and Norouzzadeh [24] used an isogeometric model to examine the size-dependent static buckling reactions of circular, elliptical, and skewed nanoplates. The shape of these nanoplates determines their buckling behavior. Eringen’s nonlocal continuum theory accounts for the nonlocal effects. In Gurtin-surface Murdoch’s elasticity hypothesis, two thin nanoplate surface layers are employed to account for the surface energy effects. This accounts for the surface energy effects. Fan et al. [25] predicted the shear buckling behavior of functionally graded skew nanoplates with surface stress. In an oblique coordinate system, the higher-order shear deformation plate theory is applied to Gurtin–Murdoch surface theory of elasticity. Several homogenization approaches have been used to estimate the mechanical properties of FGM skew nanoplates, including Reuss, Voigt, Mori–Tanaka, and Hashin–Shtrikman limit models. Norouzzadeh and Ansari [26] investigated size-dependent surface stress and nonlocal influences on square and circular nanoplates. Nanoplates can be composed of FGMs with separate surface and bulk phases. Non-local and surface effects were captured by Eringen and Gurtin–Murdoch. Mori–Tanaka distribution is used to create nanoplates. A novel matrix–vector version of the governing differential equations of motion is provided. This form uses the FEA and IGA. Using the modified coupled stress theory of elasticity, Fan et al. [27] studied the porosity-dependent nonlinear large-amplitude oscillation responses of rectangular microplates with and without a central square cutout the modified couple stress theory of elasticity (MCSTE). Based on the third-order shear deformation plate model (TSDPM), size-dependent modified couple stress-based differential motion equations are derived. The porous functionally graded material (PFGM) microplate mechanical parameters are derived using a novel power-law function that combines the material gradient and porosity dependency. The NURBS-based isogeometric technique is then utilized to satisfy C¹ continuity. Microplates with rising and decreasing porosities from the top and bottom surfaces to the center exhibit the lowest and largest pair stress size dependencies. Based on the MCSTE, Qiu et al. [28] conducted a porosity-dependent nonlinear postbuckling study on microplates constructed from PFGM. Modified coupling stress-based nonlinear differential equations were generated using the TSDPM. The effective mechanical properties of PFGM microplates can be derived using a power-law function that considers the porosity and material gradient. Next, NURBS discretization was applied to fulfill C¹ continuity criteria. Furthermore, the post-buckling phase reduced the gap between equilibrium paths for different pair stress length scales. Pham et al. [29] presented the first IGA of BDFG rectangular plates in a fluid medium. The thickness and length of the BDFG plate materials fluctuated according to the Mori–Tanaka model and power-law distributions. Rahmouni et al. [30] assessed the stress distribution and first-ply failure strength of hybrid laminated composites subjected to uniaxial tensile and compression stresses using IGA based on NURBS technology. More details on IGA and NURBS technology are available in recent studies [31–39].

The IGA method is undoubtedly one of the most reliable approaches for estimating the mechanical behavior of a structure. This method has been shown to be successful over the years, and the field of mechanics continues to show interest in further developments and improvements. In tandem with advancements in computer science, a new materials technology sector has been established. Numerous contemporary industries are progressively focusing their research and development efforts on constructing micro- and nanoscale structures. In the study, Biot’s linear poroelasticity theory with nonlocal theory was used to evaluate static bending and free and forced vibrations of functionally graded fluid-infiltrated porous (FGFP) nanoplates resting on an elastic foundation. To the best of our knowledge, there have been no published studies on this combination. With respect to the methodology, the author used the four-variable shear strain theory and isogeometric approach for skew-nanoplate and elliptical nanoplate geometries. The four-variable shear strain theory was selected over the two-variable improvement theory because it allows for complete displacements of the plate’s length, breadth, and thickness in all three dimensions. This study provides a comprehensive formula for plate calculations. In future studies, the mathematical formulae and IGA method of this study will be used to directly solve the problems of two-curvature shell structures.

The remainder of this paper is organized as follows. In Section 2, the mathematical formulations are presented, where FGFP skew-nanoplates and materials, the four-unknown shear deformation plate theory, the nonlocal constitutive relations, and the governing equation are described in detail. The isogeometric formulation of mechanical problems is introduced in Section 3. In Section 4, verification examples, numerical results, and a discussion are presented. The novel points and distributions of this study are presented in Section 5.

2 Mathematical formulation

2.1 Functionally graded fluid-infiltrated porus skew-nanoplate and materials

Fig.1 depicts the FGFP skew-nanoplate with its associated coordinate system, geometric parameters, and loading conditions. If the deviation angle is zero, then the skew nanoplate is transformed into a rectangular nanoplate. According to the power-law index, the material quality changes smoothly from the top surface (\mathit{\text{z}}=+ h/ 2) to the bottom surface (\mathit{\text{z}}=-h/2). The elastic medium is based on Pasternak’s model, which has two parameters: spring stiffness coefficient k_{\varphi } and shear stiffness coefficient k_{\mathrm{s}}.

The two master coordinate systems, namely Oxyz and the associated one O\overline{x}\overline{y}\mathit{\text{z}}, are connected to the skew plate. When Ox and O\overline{x} axes are parallel, Oy and O\overline{y} axes are at an angle \phi with respect to each other. The deviation angle illustrates the relationship between the coordinates as follows:

(1)x=\overline{x}+y\mathrm{t}\mathrm{a}\mathrm{n} \phi , y=\overline{y}/\mathrm{s}\mathrm{e}\mathrm{c} \phi .

From Eq. (1), we can determine the relationship between the first and second derivatives as follows:

\begin{array}{}\text{(2a)}& \left\{\begin{array}{c}{\displaystyle \frac{\mathrm{\partial }}{\mathrm{\partial }x}}\\ {\displaystyle \frac{\mathrm{\partial }}{\mathrm{\partial }y}}\end{array}\right\}=\left[\begin{array}{cc}1& 0\\ -\mathrm{t}\mathrm{a}\mathrm{n}\phi & \mathrm{s}\mathrm{e}\mathrm{c}\phi \end{array}\right]\left\{\begin{array}{c}{\displaystyle \frac{\mathrm{\partial }}{\mathrm{\partial }\overline{x}}}\\ {\displaystyle \frac{\mathrm{\partial }}{\mathrm{\partial }\overline{y}}}\end{array}\right\},\end{array}

\begin{array}{}\text{(2b)}& \left\{\begin{array}{c}\begin{array}{c}{\displaystyle \frac{{\mathrm{\partial }}^2}{\mathrm{\partial }{x}^2}}\\ {\displaystyle \frac{{\mathrm{\partial }}^2}{\mathrm{\partial }{y}^2}}\end{array}\\ {\displaystyle \frac{{\mathrm{\partial }}^2}{\mathrm{\partial }x\mathrm{\partial }y}}\end{array}\right\}=\left[\begin{array}{ccc}1& 0& 0\\ {\mathrm{t}\mathrm{a}\mathrm{n}}^2\phi & {\mathrm{s}\mathrm{e}\mathrm{c}}^2\phi & -2\mathrm{t}\mathrm{a}\mathrm{n}\phi \cdot \mathrm{s}\mathrm{e}\mathrm{c}\phi \\ -\mathrm{t}\mathrm{a}\mathrm{n}\phi & 0& \mathrm{s}\mathrm{e}\mathrm{c}\phi \end{array}\right]\left\{\begin{array}{c}\begin{array}{c}{\displaystyle \frac{{\mathrm{\partial }}^2}{\mathrm{\partial }{\overline{x}}^2}}\\ {\displaystyle \frac{{\mathrm{\partial }}^2}{\mathrm{\partial }{\overline{y}}^2}}\end{array}\\ {\displaystyle \frac{{\mathrm{\partial }}^2}{\mathrm{\partial }\overline{x}\mathrm{\partial }\overline{y}}}\end{array}\right\}.\end{array}

It is anticipated that the following power law can accurately depict a typical material attribute of the FGFP as a function of the plate thickness:

(3)\begin{array}{}\text{(3)}& \{\begin{array}{rl}& E\left(\mathit{\text{z}}\right)=\left[E_{\mathrm{m}}+\left(E_{\mathrm{c}}-E_{\mathrm{m}}\right){\left(\frac{1}{2}+\frac{\mathit{\text{z}}}{h}\right)}^{n_{\mathit{\text{z}}}}\right]\left(1-\lambda {\mathrm{\Gamma }}_{\mathit{\text{z}}}\right),\\ & \rho \left(\mathit{\text{z}}\right)=\left[{\rho }_{\mathrm{m}}+\left({\rho }_{\mathrm{c}}-{\rho }_{\mathrm{m}}\right){\left(\frac{1}{2}+\frac{\mathit{\text{z}}}{h}\right)}^{n_{\mathit{\text{z}}}}\right]\left(1-(1-\sqrt{1-\lambda }{)\mathrm{\Gamma }}_{\mathit{\text{z}}}\right),\end{array}\end{array}

where symbols ‘m’ and ‘c’ denote the metal and ceramic elements, respectively, \lambda denotes the coefficient of plate porosity in the range , and {\mathrm{\Gamma }}_{\mathit{\text{z}}} denotes a porosity distribution function. In this study, three different types of porosity distributions are considered [40].

Symmetric distribution (FGFP I):

\begin{array}{}\text{(4a)}& {\mathrm{\Gamma }}_{\mathit{\text{z}}}=\mathrm{c}\mathrm{o}\mathrm{s}\left(\frac{\text{π}\mathit{\text{z}}}{h}\right),\end{array}

upper asymmetric distribution (FGFP II):

\begin{array}{}\text{(4b)}& {\mathrm{\Gamma }}_{\mathit{\text{z}}}=\mathrm{c}\mathrm{o}\mathrm{s}\left(\frac{\text{π}\mathit{\text{z}}}{2h}+\frac{\text{π}}{4}\right),\end{array}

below asymmetric distribution (FGFP III):

\begin{array}{}\text{(4c)}& {\mathrm{\Gamma }}_{\mathit{\text{z}}}=\mathrm{c}\mathrm{o}\mathrm{s}\left(\frac{\text{π}\mathit{\text{z}}}{2h}-\frac{\text{π}}{4}\right).\end{array}

Given \lambda =0.5, the shape of {\mathrm{\Gamma }}_{\mathit{\text{z}}} and elastic modulus E\left(\mathit{\text{z}}\right) with respect to the power-law index n_{\mathit{\text{z}}} for three distinct types of porosity along the thickness plate is shown in Fig.2. Specifically, Fig.2(a) depicts the fluctuation of {\mathrm{\Gamma }}_{\mathit{\text{z}}} throughout the thickness of the plate and Fig.2(b)–Fig.2(d) depict the change in elastic modulus E\left(\mathit{\text{z}}\right) with respect to varying volume exponents of {\displaystyle \frac{\mathrm{A}\mathrm{l}}{{\mathrm{A}\mathrm{l}}_2{\mathrm{O}}_3}} material.

2.2 Four-unknown shear deformation plate theory

Based on the four-unknown shear deformation plate theory (FUSDT), proposed by Ref. [41], the displacement field of FGFP nanoplates is examined in relation to their thicknesses as follows:

(5)\begin{array}{}\text{(5)}& \{\begin{array}{rl}& u_x\left(x,y,\mathit{\text{z}},t\right)=u_0\left(x,y,t\right)-\mathit{\text{z}}{\varphi }_{,x}^{\mathrm{b}}-f\left(\mathit{\text{z}}\right){\varphi }_{,x}^{\mathrm{s}},\\ & u_y\left(x,y,\mathit{\text{z}},t\right)=v_0\left(x,y,t\right)-\mathit{\text{z}}{\varphi }_{,y}^{\mathrm{b}}-f\left(\mathit{\text{z}}\right){\varphi }_{,y}^{\mathrm{s}},\\ & u_{\mathit{\text{z}}}\left(x,y,\mathit{\text{z}},t\right)={\varphi }^{\mathrm{b}}+{\varphi }^{\mathrm{s}},\end{array}\end{array}

where

(6)\begin{array}{}\text{(6)}& f\left(\mathit{\text{z}}\right)=\mathit{\text{z}}-\psi \left(\mathit{\text{z}}\right),\psi \left(\mathit{\text{z}}\right)=\frac{h}{\text{π}}\mathrm{s}\mathrm{i}\mathrm{n}\frac{\text{π}\mathit{\text{z}}}{h},\end{array}

u_0 and v_0 denote the midplane displacement components along x- and y-axis, respectively, and {\varphi }^{\mathrm{b}} and {\varphi }^{\mathrm{s}} denote the bending and shear components of the transverse displacement, respectively.

The strain components of an FGFP nanoplate are obtained as follows:

\begin{array}{}\text{(7a)}& {\epsilon }_{xx}=u_{0,x}-\mathit{\text{z}}{\varphi }_{,xx}^{\mathrm{b}}-f\left(\mathit{\text{z}}\right){\varphi }_{,xx}^{\mathrm{s}},\end{array}

\begin{array}{}\text{(7b)}& {\epsilon }_{yy}=v_{0,y}-\mathit{\text{z}}{\varphi }_{,yy}^{\mathrm{b}}-f\left(\mathit{\text{z}}\right){\varphi }_{,yy}^{\mathrm{s}},\end{array}

\begin{array}{}\text{(7c)}& {\epsilon }_{xy}=u_{0,y}+v_{0,x}-2\mathit{\text{z}}{\varphi }_{,xy}^{\mathrm{b}}-2f\left(\mathit{\text{z}}\right){\varphi }_{,xy}^{\mathrm{s}},\end{array}

\begin{array}{}\text{(7d)}& {\gamma }_{x\mathit{\text{z}}}=g\left(\mathit{\text{z}}\right){\varphi }_{,x}^{\mathrm{s}},\end{array}

\begin{array}{}\text{(7e)}& {\gamma }_{y\mathit{\text{z}}}=g\left(\mathit{\text{z}}\right){\varphi }_{,y}^{\mathrm{s}},\end{array}

where

(8)\begin{array}{}\text{(8)}& g\left(\mathit{\text{z}}\right)=\mathrm{c}\mathrm{o}\mathrm{s}\frac{\text{π}\mathit{\text{z}}}{h}.\end{array}

These formulae, when written in vector form, can be shortened as follows:

(9)\begin{array}{}\text{(9)}& \begin{array}{rl}& {\mathit{\epsilon }}^{\mathrm{p}}=\left\{\begin{array}{c}{\epsilon }_{xx}\\ {\epsilon }_{yy}\\ {\epsilon }_{xy}\end{array}\right\}=\left[\begin{array}{ccc}1& \mathit{\text{z}}& f\left(\mathit{\text{z}}\right)\end{array}\right]\left\{\begin{array}{c}{\mathit{\epsilon }}^0\\ {\mathit{\epsilon }}^1\\ {\mathit{\epsilon }}^2\end{array}\right\}=\left[\begin{array}{ccc}1& \mathit{\text{z}}& f\left(\mathit{\text{z}}\right)\end{array}\right]\mathit{\epsilon },\\ & \mathit{\gamma }=\left\{\begin{array}{c}{\gamma }_{x\mathit{\text{z}}}\\ {\gamma }_{y\mathit{\text{z}}}\end{array}\right\}=g\left(\mathit{\text{z}}\right){\mathit{\gamma }}^0,\end{array}\end{array}

where

(10)\begin{array}{}\text{(10)}& \begin{array}{rl}& {\mathit{\epsilon }}^{\mathrm{p}}=\left\{\begin{array}{c}{\mathit{\epsilon }}^0\\ {\mathit{\epsilon }}^1\\ {\mathit{\epsilon }}^2\end{array}\right\},{\mathit{\epsilon }}^0=\left\{\begin{array}{c}{u}_{0,x}\\ v_{0,y}\\ u_{0,y}+v_{0,x}\end{array}\right\},{\mathit{\epsilon }}^1=-\left\{\begin{array}{c}{\varphi }_{,xx}^{\mathrm{b}}\\ {\varphi }_{,yy}^{\mathrm{b}}\\ {2\varphi }_{,xy}^{\mathrm{b}}\end{array}\right\},\\ & {\mathit{\epsilon }}^2=-\left\{\begin{array}{c}{\varphi }_{,xx}^{\mathrm{s}}\\ {\varphi }_{,yy}^{\mathrm{s}}\\ {2\varphi }_{,xy}^{\mathrm{s}}\end{array}\right\},{\mathit{\gamma }}^0=\left\{\begin{array}{c}{\varphi }_{,x}^{\mathrm{s}}\\ {\varphi }_{,y}^{\mathrm{s}}\end{array}\right\}.\end{array}\end{array}

2.3 Nonlocal constitutive relationships

According to the nonlocal theory [41,42] and linear poroelasticity theory of Biot [43], which has two features, when the porous pressure increases, the porosity increases, and when the porosity decreases, the porous pressure increases. The constitutive relations can be expressed as follows:

(11)\begin{array}{}\text{(11)}& {\sigma }_{ij}-{\mu }^2{\mathrm{\nabla }}^2{\sigma }_{ij}=2G{\epsilon }_{ij}+{\displaystyle \frac{2G{v}_{\mathrm{u}}}{1-v_{\mathrm{u}}}}\epsilon {\delta }_{ij}-\tilde{p}\alpha {\delta }_{ij},\end{array}

(12)\begin{array}{}\text{(12)}& \tilde{p}=\tilde{M}\left(\vartheta -\alpha \epsilon \right),\end{array}

(13)\begin{array}{}\text{(13)}& \tilde{M}={\displaystyle \frac{2G\left(v_{\mathrm{u}}-v\right)}{{\alpha }^2\left(1-2v_{\mathrm{u}}\right)\left(1-2v\right)}},\end{array}

(14)\begin{array}{}\text{(14)}& v_{\mathrm{u}}={\displaystyle \frac{v+{\displaystyle \frac{\alpha \beta \left(1-2v\right)}{3}}}{1-{\displaystyle \frac{\alpha \beta \left(1-2v\right)}{3}}}},\end{array}

where \tilde{p} denotes the pore fluid pressure, G denotes the shear modulus, and {\sigma }_{ij} and {\epsilon }_{ij} denote the stress and strain components, respectively. Parameter \tilde{M} denotes the Biot modulus, which is defined as the increase in the amount of fluid, v_{\mathrm{u}} denotes the undrained Poisson’s ratio , and ε denotes the volumetric strain, {\delta }_{ij} denotes the Kronecker delta function, \vartheta denotes the variation in the fluid volume content inside the pores, \alpha denotes the Biot coefficient of effective stress , β denotes the Skempton coefficient, \mu (nm) denotes the small scale effect in nanostructures, and {\mathrm{\nabla }}^2={\mathrm{\partial }}_{,xx}^2+{\mathrm{\partial }}_{,yy}^2 is the Laplacian operator. The constitutive relations of the nonlocal elasticity-stress nanoplates are finally obtained by substituting Eq. (10) in Eq. (11).

(15)\begin{array}{}\text{(15)}& \begin{array}{rl}\left(1-{\mu }^2{\mathrm{\nabla }}^2\right)& \left\{\begin{array}{c}{\sigma }_{xx}\\ {\sigma }_{yy}\\ \begin{array}{c}\begin{array}{c}{\sigma }_{xy}\\ {\sigma }_{x\mathit{\text{z}}}\end{array}\\ {\sigma }_{y\mathit{\text{z}}}\end{array}\end{array}\right\}=\left[\begin{array}{ccccc}{C}_{11}& C_{12}& 0& 0& 0\\ C_{12}& C_{22}& 0& 0& 0\\ 0& 0& C_{66}& 0& 0\\ 0& 0& 0& C_{55}& 0\\ 0& 0& 0& 0& C_{44}\end{array}\right]\\ & \cdot \left\{\begin{array}{c}\begin{array}{c}{\epsilon }_{xx}\\ {\epsilon }_{yy}\\ {\epsilon }_{xy}\end{array}\\ {\epsilon }_{x\mathit{\text{z}}}\\ {\epsilon }_{y\mathit{\text{z}}}\end{array}\right\}=\left[\begin{array}{cc}{\mathit{Q}}_{\mathrm{b}}& \mathbf{0}\\ \mathbf{0}& {\mathit{Q}}_{\mathrm{s}}\end{array}\right]\left\{\begin{array}{c}\mathit{\epsilon }\\ \mathit{\gamma }\end{array}\right\},\end{array}\end{array}

where C_{mn} is as follows:

\begin{array}{}\text{(16a)}& {\mathit{Q}}_{\mathrm{b}}=\left[\begin{array}{ccc}{C}_{11}& C_{12}& 0\\ C_{21}& C_{22}& 0\\ 0& 0& C_{66}\end{array}\right],{\mathit{Q}}_{\mathrm{s}}=\left[\begin{array}{cc}{C}_{55}& 0\\ 0& C_{44}\end{array}\right],\end{array}

\begin{array}{rl}{C}_{11}=& C_{22}={\displaystyle \frac{E\left(\mathit{\text{z}}\right)}{2\left(1+v\right)}}\frac{2}{1-v_{\mathrm{u}}^2}\\ & \cdot \left(1+v_{\mathrm{u}}+{\displaystyle \frac{\left(v_{\mathrm{u}}-v\right)\left(1+v_{\mathrm{u}}\right)}{1-2v}}\left(1-{\displaystyle \frac{S_2}{S_1}}\right)\right),\end{array}

C_{12}= \frac{E\left(\mathit{\text{z}}\right)}{2\left(1+v\right)}\frac{2}{1-v_{\mathrm{u}}^2}\left(\left(1+v_{\mathrm{u}}\right)v_{\mathrm{u}}+\frac{\left(v_{\mathrm{u}}-v\right)\left(1+v_{\mathrm{u}}\right)}{1-2v}\left(1-\frac{S_2}{S_1}\right)\right),

\begin{array}{}\text{(16b)}& C_{66}=C_{55}=C_{44}=\frac{E\left(\mathit{\text{z}}\right)}{2\left(1+v\right)},\end{array}

S_1= \frac{E\left(\mathit{\text{z}}\right)}{1+v}\left(1+\frac{v_{\mathrm{u}}}{1-2v_{\mathrm{u}}}+\frac{\left(v_{\mathrm{u}}-v\right)\left(1+v_{\mathrm{u}}\right)}{1-2v}\right),

\begin{array}{}\text{(16c)}& S_2=\frac{E\left(\mathit{\text{z}}\right)}{1+v}\left(\frac{v_{\mathrm{u}}}{1-2v_{\mathrm{u}}}+\frac{\left(v_{\mathrm{u}}-v\right)\left(1+v_{\mathrm{u}}\right)}{1-2v}\right).\end{array}

2.4 Governing equation

Using Hamilton’s principle, the general motion of the equations of the FGFP nanoplate resting on EF has been determined [43,45,46]:

(17)\begin{array}{}\text{(17)}& {\int }_{t_1}^{t_2}\left(\delta U+\delta V+\delta W-\delta T\right)\mathrm{d}t=0,\end{array}

where \delta U, \delta V, \delta W, and \delta T denote the variation components defined as follows.

The variation in the strain energy \delta U is computed as follows:

(18)\begin{array}{}\text{(18)}& \begin{array}{rl}\delta U=& {\displaystyle \underset{{\mathrm{\Omega }}_{\mathrm{e}}}{\int }{{\displaystyle \int }}_{-\frac{h}{2}}^{\frac{h}{2}}({\sigma }_{xx}\delta {\epsilon }_{xx}+{\sigma }_{yy}\delta {ϵ}_{yy}+{\sigma }_{xy}\delta {\epsilon }_{xy}+{\sigma }_{x\mathit{\text{z}}}\delta {\gamma }_{x\mathit{\text{z}}}}\\ & +{\sigma }_{y\mathit{\text{z}}}\delta {\gamma }_{y\mathit{\text{z}}})\mathrm{d}\mathit{\text{z}}\mathrm{d}x\mathrm{d}y\\ =& \underset{{\mathrm{\Omega }}_{\mathrm{e}}}{\stackrel{}{\int }}({\mathit{N}}_{xx}\delta u_{0,x}-{\mathit{M}}_{xx}\delta {\varphi }_{,xx}^{\mathrm{b}}-{\mathit{L}}_{xx}\delta {\varphi }_{,xx}^{\mathrm{s}}+{\mathit{N}}_{yy}\delta v_{0,y}\\ & -{\mathit{M}}_{yy}\delta {\varphi }_{,yy}^{\mathrm{b}}-{\mathit{L}}_{yy}\delta {\varphi }_{,yy}^{\mathrm{s}}+{\mathit{N}}_{xy}\left(\delta u_{0,y}+\delta v_{0,x}\right)-{2\mathit{M}}_{xy}\delta {\varphi }_{,xy}^{\mathrm{b}}\\ & -{2\mathit{L}}_{xy}\delta {\varphi }_{,xy}^{\mathrm{s}}+{\mathit{Q}}_{x\mathit{\text{z}}}\delta {\varphi }_{,x}^{\mathrm{s}}+{\mathit{Q}}_{y\mathit{\text{z}}}\delta {\varphi }_{,y}^{\mathrm{s}})\mathrm{d}x\mathrm{d}y.\end{array}\end{array}

The area of the element is denoted by {\mathrm{\Omega }}_{\mathrm{e}}, and the stress resultants are calculated as follows:

(19)\begin{array}{}\text{(19)}& \left\{{\mathit{N}}_{ij};{\mathit{M}}_{ij};{\mathit{L}}_{ij}\right\}={{\displaystyle \int }}_{-\frac{h}{2}}^{\frac{h}{2}}{\sigma }_{ij}\left\{1;\overline{\mathit{\text{z}}};\overline{f}\right\}\mathrm{d}\mathit{\text{z}},ij=xx,yy,xy,\end{array}

(20)\begin{array}{}\text{(20)}& {\mathit{Q}}_{ij}={{\displaystyle \int }}_{-\frac{h}{2}}^{\frac{h}{2}}{\sigma }_{ij}\overline{g}\mathrm{d}\mathit{\text{z}},ij=x\mathit{\text{z}},y\mathit{\text{z}},\end{array}

where \overline{\mathit{\text{z}}}=\mathit{\text{z}}-t_0, \overline{f}=f\left(\overline{\mathit{\text{z}}}\right), \overline{g}=g\left(\overline{\mathit{\text{z}}}\right), and t_0 denotes the distance from the mean plane to the neutral plane of the FGFP nanoplate. The distance is defined as follows:

(21)t_0= {\displaystyle \frac{{{\displaystyle \int }}_{-\frac{h}{2}}^{\frac{h}{2}}E\left(\mathit{\text{z}}\right)\mathit{\text{z}}\mathrm{d}\mathit{\text{z}}}{{{\displaystyle \int }}_{-\frac{h}{2}}^{\frac{h}{2}}E\left(\mathit{\text{z}}\right)\mathrm{d}\mathit{\text{z}}}}.

The stress resultants in terms of displacement fields (u_0,v_0,{\varphi }^{\mathrm{b}} and {\varphi }^{\mathrm{s}}) for a nonlocal model can be introduced when Eq. (15) is substituted into Eqs. (19) and (20), respectively, as follows:

(22){\left\{{\mathit{N}}_{xx}{\mathit{N}}_{yy}{\mathit{N}}_{xy}\right\}}^{\mathrm{T}}-{\mu }^2{\mathrm{\nabla }}^2{\left\{{\mathit{N}}_{xx}{\mathit{N}}_{yy}{\mathit{N}}_{xy}\right\}}^{\mathrm{T}}= {\mathit{A}\mathit{\epsilon }}^0+{\mathit{B}}^{\mathrm{b}}{\mathit{\epsilon }}^1+{\mathit{B}}^{\mathrm{s}}{\mathit{\epsilon }}^2,

(23){\left\{{\mathit{M}}_{xx}{\mathit{M}}_{yy}{\mathit{M}}_{xy}\right\}}^{\mathrm{T}}-{\mu }^2{\mathrm{\nabla }}^2{\left\{{\mathit{M}}_{xx}{\mathit{M}}_{yy}{\mathit{M}}_{xy}\right\}}^{\mathrm{T}}={{\mathit{B}}^{\mathrm{b}}\mathit{\epsilon }}^0+{\mathit{F}}^{\mathrm{b}}{\mathit{\epsilon }}^1+{\mathit{F}}^{\mathrm{s}}{\mathit{\epsilon }}^2,

(24){\left\{{\mathit{L}}_{xx}{\mathit{L}}_{yy}{\mathit{L}}_{xy}\right\}}^{\mathrm{T}}-{\mu }^2{\mathrm{\nabla }}^2{\left\{{\mathit{L}}_{xx}{\mathit{L}}_{yy}{\mathit{L}}_{xy}\right\}}^{\mathrm{T}}={{\mathit{B}}^{\mathrm{s}}\mathit{\epsilon }}^0+{\mathit{B}}^{\mathrm{s}}{\mathit{\epsilon }}^1+\mathit{H}{\mathit{\epsilon }}^2,

(25){\left\{{\mathit{Q}}_{x\mathit{\text{z}}}{\mathit{Q}}_{y\mathit{\text{z}}}\right\}}^{\mathrm{T}}-{{\mu }^2\mathrm{\nabla }}^2{\left\{{\mathit{Q}}_{x\mathit{\text{z}}}{\mathit{Q}}_{y\mathit{\text{z}}}\right\}}^{\mathrm{T}}={\mathit{A}}^{\mathrm{s}}{\mathit{\gamma }}^0,

where

(26)\begin{array}{rl}\left(\mathit{A},{\mathit{B}}^{\mathrm{b}},{\mathit{B}}^{\mathrm{s}},{\mathit{F}}^{\mathrm{b}},{\mathit{F}}^{\mathrm{s}},\mathit{H}\right)=& {\int }_{\frac{-h}{2}}^{\frac{h}{2}}{\mathit{Q}}_{\mathrm{b}}(1,\overline{\mathit{\text{z}}},\overline{f},{\overline{\mathit{\text{z}}}}^2,\overline{\mathit{\text{z}}}\overline{f},{\overline{f}}^2)\mathrm{d}\mathit{\text{z}},\\ {\mathit{A}}^{\mathrm{s}}=& {\int }_{\frac{-h}{2}}^{\frac{h}{2}}{\mathit{Q}}_{\mathrm{s}}{\overline{g}}^2\mathrm{d}\mathit{\text{z}}.\end{array}

The change in the amount of energy that is distorted in the elastic foundation \delta V is as follows:

(27)\delta V=\underset{{\mathrm{\Omega }}_{\mathrm{e}}}{\stackrel{}{\int }} \left(k_{\varphi }\left({\varphi }^{\mathrm{b}}+{\varphi }^{\mathrm{s}}\right)-k_{\mathrm{s}}{\mathrm{\nabla }}^2\left({\varphi }^{\mathrm{b}}+{\varphi }^{\mathrm{s}}\right)\right)\delta \left({\varphi }^{\mathrm{b}}+{\varphi }^{\mathrm{s}}\right)\mathrm{d} x\mathrm{d}y.

The coefficient of variation of the work done by the applied force \delta W is calculated as follows:

(28)\delta W=- \underset{{\mathrm{\Omega }}_{\mathrm{e}}}{\stackrel{}{\int }} q\left(x,y\right)\delta \left({\varphi }^{\mathrm{b}}+{\varphi }^{\mathrm{s}}\right)\mathrm{d} x\mathrm{d}y,

where q \left(x,y\right) denotes the transverse forces.

Additionally, the change in kinetic energy, denoted by \delta T, can be calculated as follows:

(29)\begin{array}{rl}\delta T=& \underset{{\mathrm{\Omega }}_{\mathrm{e}}}{\stackrel{}{\int }}{\int }_{-\frac{h}{2}}^{\frac{h}{2}}\rho \left(\mathit{\text{z}}\right)\left({\dot{u}}_x\delta {\dot{u}}_x+{\dot{u}}_y\delta {\dot{u}}_y+{\dot{u}}_{\mathit{\text{z}}}\delta {\dot{u}}_{\mathit{\text{z}}}\right)\mathrm{d}\mathit{\text{z}}\mathrm{d}x\mathrm{d}y\\ =& {\int }_{{\mathrm{\Omega }}_{\mathrm{e}}}^{}(I_0\left({\dot{u}}_0\delta {\dot{u}}_0+{\dot{v}}_0\delta {\dot{v}}_0+\left({\dot{\varphi }}^{\mathrm{b}}+{\dot{\varphi }}^{\mathrm{s}}\right)\delta \left({\dot{\varphi }}^{\mathrm{b}}+{\dot{\varphi }}^{\mathrm{s}}\right)\right)\\ & -I_1\left({\dot{u}}_0\delta {\dot{\varphi }}_{,x}^{\mathrm{b}}+{\dot{\varphi }}_{,x}^{\mathrm{b}}\delta {\dot{u}}_0+{\dot{v}}_0\delta {\dot{\varphi }}_{,y}^{\mathrm{b}}+{\dot{\varphi }}_{,y}^{\mathrm{b}}\delta {\dot{v}}_0\right)+I_2({\dot{\varphi }}_{,x}^{\mathrm{b}}\delta {\dot{\varphi }}_{,x}^{\mathrm{b}}\\ & +{\dot{\varphi }}_{,y}^{\mathrm{b}}\delta {\dot{\varphi }}_{,y}^{\mathrm{b}})-J_1\left({\dot{u}}_0\delta {\dot{\varphi }}_{,x}^{\mathrm{s}}+{\dot{\varphi }}_{,x}^{\mathrm{s}}\delta {\dot{u}}_0+{\dot{v}}_0\delta {\dot{\varphi }}_{,y}^{\mathrm{s}}+{\dot{\varphi }}_{,y}^{\mathrm{s}}\delta {\dot{v}}_0\right)\\ & +K_2\left({\dot{\varphi }}_{,x}^{\mathrm{s}}\delta {\dot{\varphi }}_{,x}^{\mathrm{s}}+{\dot{\varphi }}_{,y}^{\mathrm{s}}\delta {\dot{\varphi }}_{,y}^{\mathrm{s}}\right)+J_2({\dot{\phi }}_{,x}^{\mathrm{b}}\delta {\dot{\phi }}_{,x}^{\mathrm{s}}+{\dot{\varphi }}_{,x}^{\mathrm{s}}\delta {\dot{\varphi }}_{,x}^{\mathrm{b}}\\ & +{\dot{\varphi }}_{,y}^{\mathrm{b}}\delta {\dot{\varphi }}_{,y}^{\mathrm{s}}+{\dot{\varphi }}_{,y}^{\mathrm{s}}\delta {\dot{\varphi }}_{,y}^{\mathrm{b}}))\mathrm{d}x\mathrm{d}y,\end{array}

where (·) denotes the derivative notation over time, and I_i,J_j, and K_2 denote the mass moments of inertia defined by the following expression:

(30){(I}_0,I_1,I_2,J_1,J_2,K_2)= {\int }_{-\frac{h}{2}}^{\frac{h}{2}}\rho \left(\mathit{\text{z}}\right)\left(1,\overline{\mathit{\text{z}}},{\overline{\mathit{\text{z}}}}^2,\overline{f},\overline{\mathit{\text{z}}}\overline{f},{\overline{f}}^2\right)\mathrm{d} \mathit{\text{z}}.

By substituting Eqs. (18) and (27)−(29) into Eq. (17) and integrating by parts and collecting the coefficients of \delta u_0, {\delta v}_0,{\delta \varphi }^{\mathrm{b}} and {\delta \varphi }^{\mathrm{s}}, as provided below, the equations of motion for the FGFP nanoplates can be deduced as follows:

(31){\delta u}_0:{\mathit{N}}_{xx,x}+{\mathit{N}}_{xy,y}=I_0{\ddot{u}}_0- I_1{\ddot{\varphi }}_{,x}^{\mathrm{b}}-J_1{\ddot{\varphi }}_{,x}^{\mathrm{s}},

(32){\delta v}_0:{\mathit{N}}_{yy,y}+{\mathit{N}}_{xy,x}=I_0{\ddot{v}}_0- I_1{\ddot{\varphi }}_{,y}^{\mathrm{b}}-J_1{\ddot{\varphi }}_{,y}^{\mathrm{s}},

(33)\begin{array}{rl}& {\delta \varphi }^{\mathrm{b}}:{\mathit{M}}_{xx,xx}+{2\mathit{M}}_{xy,xy}+{\mathit{M}}_{yy,yy}+q\left(x,y\right)-k_w\left({\varphi }^{\mathrm{b}}+{\varphi }^{\mathrm{s}}\right)\\ & +k_s{\mathrm{\nabla }}^2\left({\varphi }^{\mathrm{b}}+{\varphi }^{\mathrm{s}}\right)=I_0\left({\ddot{\varphi }}^{\mathrm{b}}+{\ddot{\varphi }}^{\mathrm{s}}\right)+I_1\left({\ddot{u}}_{0,x}{+\ddot{v}}_{0,y}\right)\\ & -I_2{\mathrm{\nabla }}^2{\ddot{\varphi }}^{\mathrm{b}}-J_2{\mathrm{\nabla }}^2{\ddot{\varphi }}^{\mathrm{s}},\end{array}

(34)\begin{array}{rl}& {\delta \varphi }^{\mathrm{s}}:{\mathit{L}}_{xx,xx}+{2\mathit{L}}_{xy,xy}+{\mathit{L}}_{yy,yy}+{\mathit{Q}}_{x\mathit{\text{z}},x}+{\mathit{Q}}_{y\mathit{\text{z}},y}+q\left(x,y\right)-k_{\mathrm{w}}\left({\varphi }^{\mathrm{b}}+{\varphi }^{\mathrm{s}}\right)\\ & +k_{\mathrm{s}}{\mathrm{\nabla }}^2\left({\varphi }^{\mathrm{b}}+{\varphi }^{\mathrm{s}}\right)=I_0\left({\ddot{\varphi }}^{\mathrm{b}}+{\ddot{\varphi }}^{\mathrm{s}}\right)+J_1\left({\ddot{u}}_{0,x}{+\ddot{v}}_{0,y}\right)\\ & -J_2{\mathrm{\nabla }}^2{\ddot{\varphi }}^{\mathrm{b}}-K_2{\mathrm{\nabla }}^2{\ddot{\varphi }}^{\mathrm{s}}.\end{array}

Following this, the equations of motion for the FGFP nanoplate can be recast by inserting Eqs. (22)–(25) into Eqs. (31)–(34) as follows:

(35){\mathit{N}}_{xx,x}+ {\mathit{N}}_{xy,y}= \left(1-{\mu }^2{\mathrm{\nabla }}^2\right)\left(I_0{\ddot{u}}_0-I_1{\ddot{\varphi }}_{,x}^{\mathrm{b}}-J_1{\ddot{\varphi }}_{,x}^{\mathrm{s}}\right),

(36){\mathit{N}}_{yy,y}+ {\mathit{N}}_{xy,x}= \left(1-{\mu }^2{\mathrm{\nabla }}^2\right)\left(I_0{\ddot{v}}_0-I_1{\ddot{\varphi }}_{,y}^{\mathrm{b}}-J_1{\ddot{\varphi }}_{,y}^{\mathrm{s}}\right),

(37)\begin{array}{rl}& {\mathit{M}}_{xx,xx}+{2\mathit{M}}_{xy,xy}+{\mathit{M}}_{yy,yy}+\left(1-{\mu }^2{\mathrm{\nabla }}^2\right)q\left(x,y\right)\\ & -k_{\varphi }\left(1-{\mu }^2{\mathrm{\nabla }}^2\right)\left({\varphi }^{\mathrm{b}}+{\varphi }^{\mathrm{s}}\right)+k_{\mathrm{s}}\left(1-{\mu }^2{\mathrm{\nabla }}^2\right){\mathrm{\nabla }}^2\left({\varphi }^{\mathrm{b}}+{\varphi }^{\mathrm{s}}\right)\\ =& \left(1-{\mu }^2{\mathrm{\nabla }}^2\right)(I_0\left({\ddot{\varphi }}^{\mathrm{b}}+{\ddot{\varphi }}^{\mathrm{s}}\right)+I_1\left({\ddot{u}}_{0,x}{+\ddot{v}}_{0,y}\right)\\ & -I_2{\mathrm{\nabla }}^2{\ddot{\varphi }}^{\mathrm{b}}-J_2{\mathrm{\nabla }}^2{\ddot{\varphi }}^{\mathrm{s}}),\end{array}

(38)\begin{array}{rl}& {\mathit{L}}_{xx,xx}+{2\mathit{L}}_{xy,xy}+{\mathit{L}}_{yy,yy}+{\mathit{Q}}_{x\mathit{\text{z}},x}+{\mathit{Q}}_{y\mathit{\text{z}},y}+\left(1-{\mu }^2{\mathrm{\nabla }}^2\right)q\left(x,y\right)\\ & -k_{\phi }\left(1-{\mu }^2{\mathrm{\nabla }}^2\right)\left({\varphi }^{\mathrm{b}}+{\varphi }^{\mathrm{s}}\right)+k_{\mathrm{s}}\left(1-{\mu }^2{\mathrm{\nabla }}^2\right){\mathrm{\nabla }}^2\left({\varphi }^{\mathrm{b}}+{\varphi }^{\mathrm{s}}\right)\\ =& \left(1-{\mu }^2{\mathrm{\nabla }}^2\right)\left(I_0\left({\ddot{\varphi }}^{\mathrm{b}}+{\ddot{\varphi }}^{\mathrm{s}}\right)+J_1\left({\ddot{u}}_{0,x}{+\ddot{v}}_{0,y}\right)-J_2{\mathrm{\nabla }}^2{\ddot{\varphi }}^{\mathrm{b}}-K_2{\mathrm{\nabla }}^2{\ddot{\varphi }}^{\mathrm{s}}\right).\end{array}

Using the Galerkin technique in Eqs. (35)–(38) with the weight functions u_0, {\delta v}_0,{\delta \varphi }^{\mathrm{b}} and {\delta \varphi }^{\mathrm{s}}, respectively, and the divergence theorem, the weak form without the Neumann boundary term is as follows [46]:

(39)\begin{array}{rl}& \underset{{\mathrm{\Omega }}_{\mathrm{e}}}{\int }({\mathit{N}}_{xx}\delta u_{0,x}+{\mathit{N}}_{xy}\delta u_{0,y}+I_0\left({\ddot{u}}_0\delta u_0+{\mu }^2\left({\ddot{u}}_{0,x}\delta u_{0,x}+{\ddot{u}}_{0,y}\delta u_{0,y}\right)\right)\\ -& I_1\left({\ddot{\varphi }}_{,x}^{\mathrm{b}}\delta u_0+{\mu }^2\left({\ddot{\varphi }}_{,xx}^{\mathrm{b}}\delta u_{0,x}+{\ddot{\varphi }}_{,xy}^{\mathrm{b}}\delta u_{0,y}\right)\right)\\ -& J_1\left({\ddot{w}}_{,x}^{\mathrm{s}}\delta u_0+{\mu }^2\left({\ddot{\varphi }}_{,xx}^{\mathrm{s}}\delta u_{0,x}+{\ddot{\varphi }}_{,xy}^{\mathrm{s}}\delta u_{0,y}\right)\right))\mathrm{d}x\mathrm{d}y=0,\end{array}

(40)\begin{array}{rl}& \underset{{\mathrm{\Omega }}_{\mathrm{e}}}{\int }({\mathit{N}}_{xy}\delta v_{0,x}+{\mathit{N}}_{yy}\delta v_{0,y}+\left(I_0({\ddot{v}}_0\delta v_0+{\mu }^2\left({\ddot{v}}_{0,x}\delta v_{0,x}+{\ddot{v}}_{0,y}\delta v_{0,y}\right)\right)\\ -& I_1\left({\ddot{\varphi }}_{,y}^{\mathrm{b}}\delta v_0+{\mu }^2\left({\ddot{\varphi }}_{,xy}^{\mathrm{b}}\delta v_{0,x}+{\ddot{\varphi }}_{,yy}^{\mathrm{b}}\delta v_{0,y}\right)\right)\\ -& J_1\left({\ddot{\varphi }}_{,x}^{\mathrm{s}}\delta v_0+{\mu }^2\left({\ddot{\varphi }}_{,xx}^{\mathrm{s}}\delta v_{0,x}+{\ddot{\varphi }}_{,xy}^{\mathrm{s}}\delta v_{0,y}\right)\right))\mathrm{d}x\mathrm{d}y=0,\end{array}

(41)\begin{array}{rl}& \underset{{\mathrm{\Omega }}_{\mathrm{e}}}{\int }(-\left({\mathit{M}}_{xx}\delta {\varphi }_{,xx}^{\mathrm{b}}+2{\mathit{M}}_{xy}\delta {\varphi }_{,xy}^{\mathrm{b}}+{\mathit{M}}_{yy}\delta {\varphi }_{,yy}^{\mathrm{b}}\right)+k_{\varphi }(\left({\varphi }^{\mathrm{b}}+{\varphi }^{\mathrm{s}}\right)\delta {\varphi }^{\mathrm{b}}\\ & +{\mu }^2\left(\left({\varphi }_{,x}^{\mathrm{b}}+{\varphi }_{,x}^{\mathrm{s}}\right)\delta {\varphi }_{,x}^{\mathrm{b}}+\left({\varphi }_{,y}^{\mathrm{b}}+{\varphi }_{,y}^{\mathrm{s}}\right)\delta {\varphi }_{,y}^{\mathrm{b}}\right))+k_{\mathrm{s}}(\left({\varphi }_{,x}^{\mathrm{b}}+{\varphi }_{,x}^{\mathrm{s}}\right)\delta {\varphi }_{,x}^{\mathrm{b}}\\ & +\left({\varphi }_{,y}^{\mathrm{b}}+{\varphi }_{,y}^{\mathrm{s}}\right)\delta {\varphi }_{,y}^{\mathrm{b}}+{\mu }^2(\left({\varphi }_{,xx}^{\mathrm{b}}+w_{,xx}^{\mathrm{s}}\right)\delta {\varphi }_{,xx}^{\mathrm{b}}+\left({\varphi }_{,yy}^{\mathrm{b}}+{\varphi }_{,yy}^{\mathrm{s}}\right)\delta {\varphi }_{,yy}^{\mathrm{b}}\\ & +2\left({\phi }_{,xy}^{\mathrm{b}}+{\varphi }_{,xy}^{\mathrm{s}}\right)\delta {\varphi }_{,xy}^{\mathrm{b}}))+I_0(\left({\ddot{\varphi }}^{\mathrm{b}}+{\ddot{\varphi }}^{\mathrm{s}}\right)\delta {\varphi }^{\mathrm{b}}+{\mu }^2(\left({\ddot{\varphi }}_{,x}^{\mathrm{b}}+{\ddot{\varphi }}_{,x}^{\mathrm{s}}\right)\delta {\varphi }_{,x}^{\mathrm{b}}\\ & +\left({\ddot{\varphi }}_{,y}^{\mathrm{b}}+{\ddot{\varphi }}_{,y}^{\mathrm{s}}\right)\delta {\varphi }_{,y}^{\mathrm{b}}))-I_1({\ddot{u}}_0\delta {\varphi }_{,x}^{\mathrm{b}}+{\ddot{v}}_0\delta {\varphi }_{,y}^{\mathrm{b}}+{\mu }^2({\ddot{u}}_{0,x}\delta {\varphi }_{,xx}^{\mathrm{b}}+{\ddot{v}}_{0,x}\delta {\varphi }_{,xy}^{\mathrm{b}}\\ & +{\ddot{u}}_{0,y}\delta {\varphi }_{,xy}^{\mathrm{b}}+{\ddot{v}}_{0,y}\delta {\varphi }_{,yy}^{\mathrm{b}}))+I_2({\ddot{\varphi }}_{,x}^{\mathrm{b}}\delta {\varphi }_{,x}^{\mathrm{b}}+{\ddot{\varphi }}_{,y}^{\mathrm{b}}\delta {\varphi }_{,y}^{\mathrm{b}}+{\mu }^2({\ddot{\varphi }}_{,xx}^{\mathrm{b}}\delta {\varphi }_{,xx}^{\mathrm{b}}\\ & +2{\ddot{\varphi }}_{,xy}^{\mathrm{b}}\delta {\varphi }_{,xy}^{\mathrm{b}}+{\ddot{\varphi }}_{,yy}^{\mathrm{b}}\delta {\varphi }_{,yy}^{\mathrm{b}}))+J_2({\ddot{\varphi }}_{,x}^{\mathrm{s}}\delta {\varphi }_{,x}^{\mathrm{b}}+{\ddot{\varphi }}_{,y}^{\mathrm{s}}\delta {\varphi }_{,y}^{\mathrm{b}}+{\mu }^2({\ddot{\varphi }}_{,xx}^{\mathrm{s}}\delta {\varphi }_{,xx}^{\mathrm{b}}\\ & +2{\ddot{\varphi }}_{,xy}^{\mathrm{s}}\delta {\varphi }_{,xy}^{\mathrm{b}}+{\ddot{\varphi }}_{,yy}^{\mathrm{s}}\delta {\varphi }_{,yy}^{\mathrm{b}}))-q\left(\delta {\varphi }^{\mathrm{b}}-{\mu }^2\delta {\varphi }_{,xx}^{\mathrm{b}}-{\mu }^2\delta {\varphi }_{,yy}^{\mathrm{b}}\right))\mathrm{d}x\mathrm{d}y=0,\end{array}

(42)\begin{array}{rl}& \underset{{\mathrm{\Omega }}_{\mathrm{e}}}{\int }(-\left({\mathit{L}}_{xx}\delta {\varphi }_{,xx}^{\mathrm{s}}+2{\mathit{L}}_{xy}\delta {\varphi }_{,xy}^{\mathrm{s}}+{\mathit{L}}_{\mathit{y}\mathit{y}}\delta {\varphi }_{,yy}^{\mathrm{s}}-{\mathit{Q}}_{xz}\delta {\varphi }_{,x}^{\mathrm{s}}-{\mathit{Q}}_{yz}\delta {\varphi }_{,y}^{\mathrm{s}}\right)\\ +& k_{\varphi }\left(\left({\varphi }^{\mathrm{b}}+{\varphi }^{\mathrm{s}}\right)\delta {\varphi }^{\mathrm{s}}+{\mu }^2\left(\left({\varphi }_{,x}^{\mathrm{b}}+{\varphi }_{,x}^{\mathrm{s}}\right)\delta {\varphi }_{,x}^{\mathrm{s}}+\left({\varphi }_{,y}^{\mathrm{b}}+{\varphi }_{,y}^{\mathrm{s}}\right)\delta {\varphi }_{,y}^{\mathrm{s}}\right)\right)\\ +& k_{\mathrm{s}}(\left({\varphi }_{,x}^{\mathrm{b}}+{\varphi }_{,x}^{\mathrm{s}}\right)\delta {\varphi }_{,x}^{\mathrm{s}}+\left({\varphi }_{,y}^{\mathrm{b}}+{\varphi }_{,y}^{\mathrm{s}}\right)\delta {\varphi }_{,y}^{\mathrm{s}}+{\mu }^2(\left({\varphi }_{,xx}^{\mathrm{b}}+{\varphi }_{,xx}^{\mathrm{s}}\right)\delta {\varphi }_{,xx}^{\mathrm{s}}\\ +& \left({\varphi }_{,yy}^{\mathrm{b}}+{\varphi }_{,yy}^{\mathrm{s}}\right)\delta {\varphi }_{,yy}^{\mathrm{s}}+2\left({\varphi }_{,xy}^{\mathrm{b}}+{\varphi }_{,xy}^{\mathrm{s}}\right)\delta {\varphi }_{,xy}^{\mathrm{s}}))+I_0(\left({\ddot{\varphi }}^{\mathrm{b}}+{\ddot{\varphi }}^{\mathrm{s}}\right)\delta {\varphi }^{\mathrm{s}}\\ +& {\mu }^2\left(\left({\ddot{\varphi }}_{,x}^{\mathrm{b}}+{\ddot{\varphi }}_{,x}^{\mathrm{s}}\right)\delta {\varphi }_{,x}^{\mathrm{s}}+\left({\ddot{\varphi }}_{,y}^{\mathrm{b}}+{\ddot{\varphi }}_{,y}^{\mathrm{s}}\right)\delta {\varphi }_{,y}^{\mathrm{s}}\right))-J_1({\ddot{u}}_0\delta {\varphi }_{,x}^{\mathrm{s}}+{\ddot{v}}_0\delta {\varphi }_{,y}^{\mathrm{s}}\\ +& {\mu }^2\left({\ddot{u}}_{0,x}\delta {\varphi }_{,xx}^{\mathrm{s}}+{\ddot{v}}_{0,x}\delta {\varphi }_{,xy}^{\mathrm{s}}+{\ddot{u}}_{0,y}\delta {\varphi }_{,xy}^{\mathrm{s}}+{\ddot{v}}_{0,y}\delta {\varphi }_{,yy}^{\mathrm{s}}\right))+J_2({\ddot{\varphi }}_{,x}^{\mathrm{b}}\delta {\varphi }_{,x}^{\mathrm{s}}\\ +& {\ddot{\varphi }}_{,y}^{\mathrm{b}}\delta {\varphi }_{,y}^{\mathrm{s}}+{\mu }^2\left({\ddot{\varphi }}_{,xx}^{\mathrm{b}}\delta {\varphi }_{,xx}^{\mathrm{s}}+2{\ddot{\varphi }}_{,xy}^{\mathrm{b}}\delta {\varphi }_{,xy}^{\mathrm{s}}+{\ddot{\varphi }}_{,yy}^{\mathrm{b}}\delta {\varphi }_{,yy}^{\mathrm{s}}\right))+K_2({\ddot{\varphi }}_{,x}^{\mathrm{s}}\delta {\varphi }_{,x}^{\mathrm{s}}\\ +& {\ddot{\varphi }}_{,y}^{\mathrm{s}}\delta {\varphi }_{,y}^{\mathrm{s}}+{\mu }^2\left({\ddot{\varphi }}_{,xx}^{\mathrm{s}}\delta {\varphi }_{,xx}^{\mathrm{s}}+2{\ddot{\varphi }}_{,xy}^{\mathrm{s}}\delta {\varphi }_{,xy}^{\mathrm{s}}+{\ddot{\varphi }}_{,yy}^{\mathrm{s}}\delta {\varphi }_{,yy}^{\mathrm{s}}\right))\\ & -q\left(\delta {\varphi }^{\mathrm{s}}-{\mu }^2\delta {\varphi }_{,xx}^{\mathrm{s}}-{\mu }^2\delta {\varphi }_{,yy}^{\mathrm{s}}\right))\mathrm{d}x\mathrm{d}y=0.\end{array}

When Eqs. (39)–(42) are added together, we arrive at the following equation:

(43)\begin{array}{rl}& \underset{{\mathrm{\Omega }}_{\mathrm{e}}}{\stackrel{}{\int }}\left({\mathit{\epsilon }}^{\mathrm{T}}{\mathit{D}}_{\mathrm{b}}\delta \mathit{\epsilon }+{\mathit{\gamma }}^{\mathrm{T}}{\mathit{D}}_{\mathrm{s}}\delta \mathit{\gamma }\right)\mathrm{d}x\mathrm{d}y+k_{\varphi }(\underset{{\mathrm{\Omega }}_{\mathrm{e}}}{\stackrel{}{\int }}{u}_{\mathit{\text{z}}}^{\mathrm{T}}\delta u_{\mathit{\text{z}}}\mathrm{d}x\mathrm{d}y\\ & +{\mu }^2\underset{{\mathrm{\Omega }}_{\mathrm{e}}}{\stackrel{}{\int }}\left(u_{\mathit{\text{z}},x}^{\mathrm{T}}\delta u_{\mathit{\text{z}},x}+u_{\mathit{\text{z}},y}^{\mathrm{T}}\delta u_{\mathit{\text{z}},y}\right)\mathrm{d}x\mathrm{d}y)+k_{\mathrm{s}}(\underset{{\mathrm{\Omega }}_{\mathrm{e}}}{\stackrel{}{\int }}\mathrm{\nabla }{u}_{\mathit{\text{z}}}^{\mathrm{T}}\mathrm{\nabla }\delta u_{\mathit{\text{z}}}\mathrm{d}x\mathrm{d}y\\ & +{\mu }^2\underset{{\mathrm{\Omega }}_{\mathrm{e}}}{\stackrel{}{\int }}\left({\mathrm{\nabla }u}_{\mathit{\text{z}},x}^{\mathrm{T}}\mathrm{\nabla }\delta u_{\mathit{\text{z}},x}+{\mathrm{\nabla }u}_{\mathit{\text{z}},y}^{\mathrm{T}}\mathrm{\nabla }\delta u_{\mathit{\text{z}},y}\right)\mathrm{d}x\mathrm{d}y)\\ =& \underset{{\mathrm{\Omega }}_{\mathrm{e}}}{\stackrel{}{\int }}q\left(\delta u_{\mathit{\text{z}}}-{\mu }^2\delta u_{\mathit{\text{z}},xx}-{\mu }^2\delta u_{\mathit{\text{z}},yy}\right)\mathrm{d}x\mathrm{d}y+{\int }_{{\mathrm{\Omega }}_{\mathrm{e}}}^{}{\mathit{u}}^{\mathrm{T}}{\mathit{H}}_{\mathrm{m}}\delta \ddot{\mathit{u}}\mathrm{d}x\mathrm{d}y\\ & +{\mu }^2{\int }_{{\mathrm{\Omega }}_{\mathrm{e}}}\left({\mathit{u}}_{,x}^{\mathrm{T}}{\mathit{H}}_{\mathrm{m}}\delta {\ddot{\mathit{u}}}_{,x}+{\mathit{u}}_{,y}^{\mathrm{T}}{\mathit{H}}_{\mathrm{m}}\delta {\ddot{\mathit{u}}}_{,y}\right)\mathrm{d}x\mathrm{d}y,\end{array}

where

(44){\mathit{D}}_{\mathbf{b}}=\left[\begin{array}{ccc}\mathit{A}& {\mathit{B}}^{\mathbf{b}}& {\mathit{B}}^{\mathbf{s}}\\ {\mathit{B}}^{\mathbf{b}}& {\mathit{F}}^{\mathbf{b}}& {\mathit{F}}^{\mathbf{s}}\\ {\mathit{B}}^{\mathbf{s}}& {\mathit{F}}^{\mathbf{s}}& \mathit{H}\end{array}\right], {\mathit{D}}_{\mathbf{s}}= {\mathit{A}}^{\mathbf{s}},

(45)\mathit{u}=\left\{\begin{array}{c}{u}_0\\ v_0\\ {\varphi }^{\mathrm{b}}\\ {\varphi }^{\mathrm{s}}\\ {\varphi }_{,x}^{\mathrm{b}}\\ {\varphi }_{,x}^{\mathrm{s}}\\ {\varphi }_{,y}^{\mathrm{b}}\\ {\varphi }_{,y}^{\mathrm{s}}\end{array}\right\},{\mathit{H}}_{\mathrm{m}}=\left[\begin{array}{cccccccc}{I}_0& 0& 0& 0& I_1& J_1& 0& 0\\ 0& I_0& 0& 0& 0& 0& I_1& J_1\\ 0& 0& I_0& 0& 0& 0& 0& 0\\ 0& 0& 0& I_0& 0& 0& 0& 0\\ I_1& 0& 0& 0& I_2& J_2& 0& 0\\ J_1& 0& 0& 0& J_2& K_2& 0& 0\\ 0& I_1& 0& 0& 0& 0& I_2& J_2\\ 0& J_1& 0& 0& 0& 0& J_2& K_2\end{array}\right].

3 Isogeometric formulation for mechanical problems

3.1 Non-uniform rational B-splines basis functions: Basis of isogeometric analysis

In a space with only one dimension, knot vector k\left(\zeta \right) is a set of numbers that do not decrease and are between zero and one. This set of numbers is represented by equation k\left(\zeta \right)=\left\{{\zeta }_1=0,\dots ,{\zeta }_i,\dots ,{\zeta }_{n+p+1}=1\right\}, where i denotes the knot index, {\zeta }_i denotes the ith knot, n denotes the number of basic functions, and p denotes the order of the polynomial. The following is a recursive definition of the ith B-spline basis function of degree p, which is represented by the notation N_{i,p}\left(\zeta \right) [4].

(46)

and

(47)N_{i,p}\left(\zeta \right)={\displaystyle \frac{\zeta -{\zeta }_i}{{\zeta }_{i+p}-{\zeta }_i}}{N}_{i,p-1}\left(\zeta \right)+{\displaystyle \frac{{\zeta }_{i+p+1}-\zeta }{{\zeta }_{i+p+1}-{\zeta }_{i+1}}}{N}_{i+1,p-1}\left(\zeta \right),\mathrm{f}\mathrm{o}\mathrm{r}p\ge 1.

The construction of two-dimensional NURBS basis functions involves multiplying two univariate B-spline basis functions [4] as follows:

(48)R_{i,j}^{p,q}\left(\zeta ,\eta \right)={\displaystyle \frac{N_{i,p}\left(\zeta \right)N_{j,q}\left(\eta \right)w_{i,j}}{\sum _{i=1}^n\sum _{j=1}^m{N}_{i,p}\left(\zeta \right)N_{j,q}\left(\eta \right)w_{i,j}}},

where w_{i,j} denotes weight; N_{i,p}\left(\zeta \right) and N_{j,q}\left(\eta \right) signify the B-spline basis functions of order p in \zeta direction and order p in \eta direction; N_{j,q}\left(\eta \right) follows the recursive formula described in Eq. (47) with the knot vector k\left(\eta \right), which is defined in the same manner as k\left(\zeta \right).

3.2 Non-uniform rational B-splines-based finite element formulation

The following equations can be used to calculate the plate displacements in the middle plane.

(49){\mathit{q}}^{\mathrm{h}}=\sum _{e=1}^{Ne}{\mathit{R}}_e{\mathit{q}}_e,

with

(50){\mathit{q}}_e=\left[u_{0e},v_{0e},{\varphi }_e^b,{\varphi }_e^s\right],

where Ne=\left(p+1\right)\left(q+1\right) denotes the number of control points per physical element, and {\mathit{R}}_e and q_e indicate the shape functions and unknown displacement vector at control point e, respectively.

The in-plane shear and normal stresses can be rewritten by substituting Eq. (49) into Eq. (10), which allows for the following expression.

(51)\begin{array}{rl}& {\mathit{\epsilon }}^0=\sum _{e=1}^{Ne}{\mathit{B}}_e^1{\mathit{q}}_e,{\mathit{\epsilon }}^1=\sum _{e=1}^{Ne}{\mathit{B}}_e^2{\mathit{q}}_e,\\ & {\mathit{\epsilon }}^2=\sum _{e=1}^{Ne}{\mathit{B}}_e^3{\mathit{q}}_e,{\mathit{\gamma }}^0=\sum _{e=1}^{Ne}{\mathit{B}}_e^{\mathrm{b}}{\mathit{q}}_e,\end{array}

where

(52)\begin{array}{rl}& {\mathit{B}}_s^1=\left[\begin{array}{cccc}{\mathit{R}}_{\theta ,x}& 0& 0& 0\\ 0& {\mathit{R}}_{\epsilon ,y}& 0& 0\\ {\mathit{R}}_{\theta ,v}& {\mathit{R}}_{\epsilon ,x}& 0& 0\end{array}\right],{\mathit{B}}_i^2=\left[\begin{array}{cccc}0& 0& -{\mathit{R}}_{\epsilon ,2x}& 0\\ 0& 0& -{\mathit{R}}_{\epsilon ,2w}& 0\\ 0& 0& -2{\mathit{R}}_{\epsilon ,2y}& 0\end{array}\right],\\ & {\mathit{B}}_e^3=\left[\begin{array}{cccc}0& 0& 0& {-\mathit{R}}_{e,xx}\\ 0& 0& 0& {-\mathit{R}}_{e,yy}\\ 0& 0& 0& {-2\mathit{R}}_{e,xy}\end{array}\right],{\mathit{B}}_e^{\mathrm{b}}=\left[\begin{array}{cccc}0& 0& 0& {\mathit{R}}_{e,x}\\ 0& 0& 0& {\mathit{R}}_{e,y}\end{array}\right].\end{array}

Substituting Eq. (49) into Eq. (45), the displacement vector \mathit{u} is expressed as follows:

(53)\mathit{u}=\sum _{e=1}^{Ne} {\mathit{N}}_e{\mathit{q}}_e,

with

{\mathit{N}}_e=\left[\begin{array}{ccc}{\mathit{N}}_e^1;& {\mathit{N}}_e^2;& {\mathit{N}}_e^3\end{array}\right],

(54)\begin{array}{rl}& {\mathit{N}}_e^1=\left[\begin{array}{cccc}{\mathit{R}}_e& 0& 0& 0\\ 0& 0& {\mathit{R}}_{e,x}& 0\\ 0& 0& 0& {\mathit{R}}_{e,x}\end{array}\right],{\mathit{N}}_e^2=\left[\begin{array}{cccc}0& {\mathit{R}}_e& 0& 0\\ 0& 0& {\mathit{R}}_{e,y}& 0\\ 0& 0& 0& {\mathit{R}}_{e,y}\end{array}\right],\\ & {\mathit{N}}_e^3=\left[\begin{array}{cccc}0& 0& {\mathit{R}}_e& 0\\ 0& 0& 0& {\mathit{R}}_e\\ 0& 0& 0& 0\end{array}\right].\end{array}

The displacement vector u_{\mathit{\text{z}}} is expressed as follows:

(55)u_{\mathit{\text{z}}}= \sum _{e=1}^{Ne}\left[\begin{array}{cccc}0& 0& {\mathit{R}}_e& {\mathit{R}}_e\end{array}\right]{\mathit{q}}_e=\sum _{e=1}^{Ne} {\mathit{B}}_e^4{\mathit{q}}_e.

By replacing Eqs. (1) and (2) into Eqs. (51), (53), (55) the definitions of the strain matrices in Oxyz coordinate system connected to the skew plate are as follows:

\begin{array}{}\text{(56a)}& {\overline{\mathit{B}}}_e^1=\left[\begin{array}{cccc}{\mathit{R}}_{e,\overline{x}}& 0& 0& 0\\ 0& {\mathit{R}}_{e,\overline{y}}\mathrm{s}\mathrm{e}\mathrm{c}\phi -{\mathit{R}}_{e,\overline{x}}\mathrm{t}\mathrm{a}\mathrm{n}\phi & 0& 0\\ {\mathit{R}}_{e,\overline{y}}\mathrm{s}\mathrm{e}\mathrm{c}\phi -{\mathit{R}}_{e,\overline{x}}\mathrm{t}\mathrm{a}\mathrm{n}\phi & {\mathit{R}}_{e,\overline{x}}& 0& 0\end{array}\right],\end{array}

\begin{array}{}\text{(56b)}& {\overline{\mathit{B}}}_e^2=\left[\begin{array}{cccc}0& 0& {-\mathit{R}}_{e,\overline{x}\overline{x}}& 0\\ 0& 0& 2{\mathit{R}}_{e,\overline{x}\overline{y}}\mathrm{t}\mathrm{a}\mathrm{n}\phi {-\mathit{R}}_{e,\overline{x}\overline{x}}{\mathrm{t}\mathrm{a}\mathrm{n}}^2\phi {-\mathit{R}}_{e,\overline{y}\overline{y}}{\mathrm{s}\mathrm{e}\mathrm{c}}^2\phi & 0\\ 0& 0& 2\left({\mathit{R}}_{e,\overline{x}\overline{x}}\mathrm{t}\mathrm{a}\mathrm{n}\phi -{\mathit{R}}_{e,\overline{x}\overline{y}}\mathrm{s}\mathrm{e}\mathrm{c}\phi \right)& 0\end{array}\right],\end{array}

\begin{array}{}\text{(56c)}& {\overline{\mathit{B}}}_e^3=\left[\begin{array}{cccc}0& 0& 0& {-R}_{e,\overline{x}\overline{x}}\\ 0& 0& 0& 2{\mathit{R}}_{e,\overline{x}\overline{y}}\mathrm{t}\mathrm{a}\mathrm{n}\phi {-\mathit{R}}_{e,\overline{x}\overline{x}}{\mathrm{t}\mathrm{a}\mathrm{n}}^2\phi {-\mathit{R}}_{e,\overline{y}\overline{y}}{\mathrm{s}\mathrm{e}\mathrm{c}}^2\phi \\ 0& 0& 0& 2\left({\mathit{R}}_{e,\overline{x}\overline{x}}\mathrm{t}\mathrm{a}\mathrm{n}\phi -{\mathit{R}}_{e,\overline{x}\overline{y}}\mathrm{s}\mathrm{e}\mathrm{c}\phi \right)\end{array}\right],\end{array}

\begin{array}{}\text{(56d)}& {\overline{\mathit{B}}}_e^{1\mathrm{b}}=\left[\begin{array}{cccc}0& 0& 0& {\mathit{R}}_{e,\overline{x}}\\ 0& 0& 0& {\mathit{R}}_{e,\overline{y}}\mathrm{s}\mathrm{e}\mathrm{c}\phi -{\mathit{R}}_{e,\overline{x}}\mathrm{t}\mathrm{a}\mathrm{n}\phi \end{array}\right].\end{array}

\begin{array}{}\text{(57a)}& {\overline{\mathit{N}}}_e=\left[\begin{array}{ccc}{\overline{\mathit{N}}}_e^1;& {\overline{\mathit{N}}}_e^2;& {\overline{\mathit{N}}}_e^3\end{array}\right],\end{array}

\begin{array}{}\text{(57b)}& {\overline{\mathit{N}}}_e^1=\left[\begin{array}{cccc}{\mathit{R}}_e& 0& 0& 0\\ 0& 0& {\mathit{R}}_{e,\overline{x}}& 0\\ 0& 0& 0& {\mathit{R}}_{e,\overline{x}}\end{array}\right],\end{array}

\begin{array}{}\text{(57c)}& {\overline{\mathit{N}}}_e^2=\left[\begin{array}{cccc}0& {\mathit{R}}_e& 0& 0\\ 0& 0& {\mathit{R}}_{e,\overline{y}}\mathrm{s}\mathrm{e}\mathrm{c}\phi -{\mathit{R}}_{e,\overline{x}}\mathrm{t}\mathrm{a}\mathrm{n}\phi & 0\\ 0& 0& 0& {\mathit{R}}_{e,\overline{y}}\mathrm{s}\mathrm{e}\mathrm{c}\phi -{\mathit{R}}_{e,\overline{x}}\mathrm{t}\mathrm{a}\mathrm{n}\phi \end{array}\right],\end{array}

\begin{array}{}\text{(57d)}& {\overline{\mathit{N}}}_e^3=\left[\begin{array}{cccc}0& 0& {\mathit{R}}_e& 0\\ 0& 0& 0& {\mathit{R}}_e\\ 0& 0& 0& 0\end{array}\right].\end{array}

{\overline{\mathit{B}}}_e^4= \left[\begin{array}{cccc}0& 0& {\mathit{R}}_e& {\mathit{R}}_e\end{array}\right].

After substituting Eqs. (56)–(58) into Eq. (43), we can obtain the global equilibrium equations for static bending, free vibration, and forced vibration using the following formulae:

For the static bending problem:

(59)\sum _{i=1}^{Tn}\left({\mathbb{K}}_{\mathrm{p}}^i+{\mathbb{K}}_{\mathrm{f}}^i\right)\mathit{q}=\sum _{i=1}^{Tn}{\overline{\mathbb{F}}}_{\mathrm{p}}^i.

For the free vibration problem: setting \mathit{q}=\overline{\mathit{q}}\mathrm{s} \mathrm{i}\mathrm{n}\left(\omega t\right), \omega denotes the natural frequency.

(60)\left(\sum _{i=1}^{Tn}\left({\mathbb{K}}_{\mathrm{p}}^i+{\mathbb{K}}_{\mathrm{f}}^i\right)-{\omega }^2\sum _{i=1}^{Tn}{\mathbb{M}}_{\mathrm{p}}^i\right)\overline{\mathit{q}}=0.

For the forced vibration problem:

(61)\sum _{i=1}^{Tn}{\mathbb{M}}_{\mathrm{p}}^i\ddot{\mathit{q}} +\mathbb{C}\dot{\mathit{q}}+\sum _{i=1}^{Tn}\left({\mathbb{K}}_{\mathrm{p}}^i+{\mathbb{K}}_{\mathrm{f}}^i\right)\mathit{q}=\sum _{i=1}^{Tn}{\overline{\mathbb{F}}}_{\mathrm{p}}^i\left(t\right) ,

where Tn denotes the total number of nanoplate elements, and the matrices are defined as follows:

The stiffness matrix is as follows:

(62){\mathbb{K}}_{\mathrm{p}}^i={\int }_{{\mathrm{\Omega }}_{\mathrm{e}}}^{}{\left\{\begin{array}{c}{\overline{\mathit{B}}}^1\\ {\overline{\mathit{B}}}^2\\ {\overline{\mathit{B}}}^3\\ {\overline{\mathit{B}}}^{\mathrm{b}}\end{array}\right\}}^{\mathrm{T}}\left[\begin{array}{cccc}\mathit{A}& {\mathit{B}}^{\mathrm{b}}& {\mathit{B}}^{\mathrm{s}}& {\mathbf{0}}_{3x2}\\ {\mathit{B}}^{\mathrm{b}}& {\mathit{F}}^{\mathrm{b}}& {\mathit{F}}^{\mathrm{s}}& {\mathbf{0}}_{3x2}\\ {\mathit{B}}^{\mathrm{s}}& {\mathit{F}}^{\mathrm{s}}& \mathit{H}& {\mathbf{0}}_{3x2}\\ {\mathbf{0}}_{2x3}& {\mathbf{0}}_{2x3}& {\mathbf{0}}_{2x3}& {\mathit{A}}^{\mathrm{s}}\end{array}\right]\left\{\begin{array}{c}{\overline{\mathit{B}}}^1\\ {\overline{\mathit{B}}}^2\\ {\overline{\mathit{B}}}^3\\ {\overline{\mathit{B}}}^{\mathrm{b}}\end{array}\right\}\mathrm{d}\overline{x}\mathrm{d} \overline{y},

(63)\begin{array}{rl}{\mathbb{K}}_{\mathrm{f}}^i=& \underset{{\mathrm{\Omega }}_{\mathrm{e}}}{\int }(k_{\mathrm{w}}\left({\overline{\mathit{B}}}^{4\mathrm{T}}{\overline{\mathit{B}}}^4+{\mu }^2\left({\overline{\mathit{B}}}_{,\overline{x}}^{4\mathrm{T}}{\overline{\mathit{B}}}_{,\overline{x}}^4+{\overline{\mathit{B}}}_{,\overline{y}}^{4\mathrm{T}}{\overline{\mathit{B}}}_{,\overline{y}}^4\right)\right)\\ & +k_{\mathrm{s}}\left({\overline{\mathit{B}}}_{,\overline{x}}^{4\mathrm{T}}{\overline{\mathit{B}}}_{,\overline{x}}^4+{\overline{\mathit{B}}}_{,\overline{y}}^{4\mathrm{T}}{\overline{\mathit{B}}}_{,\overline{y}}^4\right)+{\mu }^2{k}_{\mathrm{s}}({\overline{\mathit{B}}}_{,\overline{x}\overline{x}}^{4\mathrm{T}}{\overline{\mathit{B}}}_{,\overline{x}\overline{x}}^4+{\overline{\mathit{B}}}_{,\overline{y}\overline{y}}^{4\mathrm{T}}{\overline{\mathit{B}}}_{,\overline{y}\overline{y}}^4\\ & +{\overline{\mathit{B}}}_{,\overline{x}\overline{x}}^{4\mathrm{T}}{\overline{\mathit{B}}}_{,\overline{y}\overline{y}}^4+{\overline{\mathit{B}}}_{,\overline{y}\overline{y}}^{4\mathrm{T}}{\overline{\mathit{B}}}_{,\overline{x}\overline{x}}^4))\mathrm{d}\overline{x}\mathrm{d}\overline{y}.\end{array}

The mass matrix is as follows:

(64){\mathbb{M}}_{\mathrm{p}}^i={\int }_{{\mathrm{\Omega }}_{\mathrm{e}}}^{}\left({\overline{\mathit{N}}}^{\mathrm{T}}{\mathit{H}}_{\mathrm{m}}\overline{\mathit{N}}+{\mu }^2\left({\overline{\mathit{N}}}_{,\overline{x}}^{\mathrm{T}}{\mathit{H}}_{\mathrm{m}}{\overline{\mathit{N}}}_{,\overline{x}}+{\overline{\mathit{N}}}_{,\overline{y}}^{\mathrm{T}}{\mathit{H}}_{\mathrm{m}}{\overline{\mathit{N}}}_{,\overline{y}}\right)\right)\mathrm{d} \overline{x}\mathrm{d} \overline{y}.

The load vector is as follows:

(65){\overline{\mathbb{F}}}_{\mathrm{p}}^i=-{\int }_{{\mathrm{\Omega }}_{\mathrm{e}}}^{}\left\{q\left(x,y,t\right)\left({\overline{\mathit{B}}}^{4\mathrm{T}}-{\mu }^2\left({\overline{\mathit{B}}}_{,xx}^{4\mathrm{T}}+{\overline{\mathit{B}}}_{,yy}^{4\mathrm{T}}\right)\right)\right\}\mathrm{d} \overline{x}\mathrm{d} \overline{y}.

The damping matrix \mathbb{C} is as follows:

(66)\mathbb{C}={ϵ}_1\sum _{i=1}^{Tn}{\mathbb{M}}_{\mathrm{p}}^i+{ϵ}_2\sum _{i=1}^{Tn}\left({\mathbb{K}}_{\mathrm{p}}^i+{\mathbb{K}}_{\mathrm{f}}^i\right),

where {ϵ}_1 and {ϵ}_2 denote the damping constants determined from the damping ratio and oscillation’s natural frequency. For the convenience of calculation, we often only examine the first two frequencies, {\omega }_1 and {\omega }_2, and under the assumption that the drag ratio is constant, {\xi }_1={\xi }_2=\xi, Rayleigh constants {ϵ}_1 and {ϵ}_2 are calculated as follows:

(67){ϵ}_1= {\displaystyle \frac{2\xi }{{\omega }_1+{\omega }_2}}{\omega }_1{\omega }_2, {ϵ}_2= {\displaystyle \frac{2\xi }{{\omega }_1+{\omega }_2}}.

Additionally, the direct integration method, termed as Newmark-beta, is applied to circumvent the problem of the nanoplate exhibiting a dynamic response when it is exposed to a dynamic load. The algorithm for this approach and its diagram are shown below.

The following is a description of the direct integration algorithm developed by Newmark.

The first step involves determining the starting condition:

(68)q_0= 0,{\dot{q}}_0=0,

(69)\overline{\mathit{A}}=\mathbb{M}+\overline{\alpha }\mathrm{d} t\mathbb{C}+\overline{\beta }{\mathrm{d}t}^2\mathbb{K},

where \mathrm{d} t denotes the integral step, and the following constants denote the difference:

(70)\overline{\alpha }=0.5,\overline{\beta }=0,25.

The following formula is used to determine the acceleration at step \left(n+1\right):

(71)\begin{array}{rl}{\ddot{\mathit{q}}}_{n+1}=& \frac{1}{\overline{\mathit{A}}}(\overline{\mathbb{F}}-\mathbb{C}\left({\dot{\mathit{q}}}_n+\left(1-\overline{\alpha }\right)\mathrm{d}t{\ddot{\mathit{q}}}_n\right)\\ & -\mathbb{K}\left({\mathit{q}}_n+\mathrm{d}t{\dot{\mathit{q}}}_n+\left(\frac{1}{2}-\overline{\beta }\right){\mathrm{d}t}^2{\ddot{\mathit{q}}}_n\right)).\end{array}

To calculate the velocity at step \left(n+1\right), the following formula is applied:

(72){\dot{\mathit{q}}}_{n+1}= {\dot{\mathit{q}}}_n+\left(1-\overline{\alpha }\right)\mathrm{d} t{\ddot{\mathit{q}}}_n+ \overline{\beta }\mathrm{d} t{\ddot{\mathit{q}}}_n.

The following equation should be used to compute the displacement at step \left(n+1\right):

(73){\mathit{q}}_{n+1}={\mathit{q}}_n+\mathrm{d}t{\dot{\mathit{q}}}_n+{\mathrm{d}t}^2\left(\frac{1}{2}-\overline{\beta }\right){\ddot{\mathit{q}}}_n+ \overline{\beta }{\mathrm{d}t}^2{\ddot{\mathit{q}}}_{n+1}.

In the following numerical examples, the following boundary conditions (BCs) are considered.

- Simply supported (S)

For upper and lower edges:u_0={\varphi }^{\mathrm{b}}={\varphi }^{\mathrm{s}}={\varphi }_{,x}^{\mathrm{b}}= {\varphi }_{,x}^{\mathrm{s}}= 0.

For oblique edges:v_0={\varphi }^{\mathrm{b}}={\varphi }^{\mathrm{s}}={\varphi }_{,\stackrel{-}{y}}^{\mathrm{b}}={\varphi }_{,\stackrel{-}{y}}^{\mathrm{s}}=0.

- Clamped (C)

For upper and lower edges:u_0=v_0={\varphi }^{\mathrm{b}}={\varphi }^{\mathrm{s}}={\varphi }_{,x}^{\mathrm{b}}= {\varphi }_{,x}^{\mathrm{s}}= {\varphi }_{,y}^{\mathrm{b}}= {\varphi }_{,y}^{\mathrm{s}}= 0.

For oblique edges:u_0=v_0={\varphi }^{\mathrm{b}}={\varphi }^{\mathrm{s}}={\varphi }_{,\overline{x}}^{\mathrm{b}}={\varphi }_{,\overline{x}}^{\mathrm{s}}={\varphi }_{,\overline{y}}^{\mathrm{b}}={\varphi }_{,\overline{y}}^{\mathrm{s}}=0.

- Free (F): all degrees of freedom (DOFs) at the boundary differ from zero.

The Dirichlet BCs u_0, v_0, {\varphi }^{\mathrm{b}}, and {\varphi }^{\mathrm{s}} can be addressed in the same manner as in the conventional FEM when performing IGA. Nevertheless, derivatives {\varphi }_{,x}^{\mathrm{b}}, {\varphi }_{,x}^{\mathrm{s}}, {\varphi }_{,\overline{y}}^{\mathrm{b}}, and {\varphi }_{,\overline{y}}^{\mathrm{s}} should be handled in a specific manner [47]. Nguyen et al. [48] provided the simplest approach, which involves imposing zero values for all of the control points’ deflections that are connected to clamped BCs and for those that are close to it. Hence, as a direct result, clamped BCs are essentially compelled to satisfy a condition that requires zero slope. The accuracy of this straightforward procedure was demonstrated by the results obtained for the clamped BCs.

4 Numerical results

In this section, MATLAB is used to create a calculation program based on the formulae developed in the previous two sections (2 and 3). The accuracy and reliability of the calculation program were verified by a comparison with the published assertions derived from the model presented in the article when applied to specific scenarios. The impact of the nonlocal coefficient, porosity coefficient, power-law index, Skempton coefficient index, elastic foundation stiffness coefficient, and other parameters on static bending, free and forced vibration of the FGFP skew, and elliptical nanoplate were examined via the application of many cases. The following is a breakdown of the mechanical properties.

\mathrm{A}\mathrm{l}:E_{\mathrm{m}}=70\mathrm{G}\mathrm{P}\mathrm{a} ,{\rho }_{\mathrm{m}}= 2700\mathrm{k}\mathrm{g}/{\mathrm{m}}^3, {\nu }_{\mathrm{m}}= 0.3,

{\mathrm{A}\mathrm{l}}_2{\mathrm{O}}_3:E_{\mathrm{c}}=380\mathrm{G}\mathrm{P}\mathrm{a} ,{\rho }_{\mathrm{c}}= 3800\mathrm{k}\mathrm{g}/{\mathrm{m}}^3, {\nu }_{\mathrm{c}}= \mathrm{0.3.}

To facilitate the comparison of the findings with previous studies, the following dimensionless formulae were introduced:

\begin{array}{rl}& w_1=\frac{100E_{\mathrm{c}}{h}^3}{q_0{a}^4}{w}_{\mathrm{c}},{\mathrm{\Omega }}_{\mathrm{m}}={\omega }_{\mathrm{m}}\frac{b^2}{{\text{π}}^2}\sqrt{\frac{{\rho }_{\mathrm{c}}h}{D_{\mathrm{c}}}},\\ & {\hat{\mathrm{\Omega }}}_1={\omega }_{1}h\sqrt{\frac{{\rho }_{\mathrm{c}}}{E_{\mathrm{c}}}},W_1=\frac{100E_{\mathrm{c}}{h}^3}{q_0{a}^4}{W}_{\mathrm{c}}\left(t\right),\end{array}

K_{\mathrm{w}}= \frac{k_{\phi }{a}^4}{D_{\mathrm{c}}},K_{\mathrm{s}}= \frac{k_{\mathrm{s}}{a}^2}{D_{\mathrm{c}}},D_{\mathrm{c}}= \frac{E_{\mathrm{c}}{h}^3}{12\left(1-{\nu }^2\right)} ,

where w_{\mathrm{c}} denotes the deflection at the nanoplate center under a uniformly distributed static load. Furthermore, W_{\mathrm{c}}\left(t\right) denotes the deflection at the skew and elliptical nanoplate center under a dynamic load.

4.1 Verification examples

In this section, the performance characteristics are evaluated, including convergence analysis and numerical validation. Tab.1 presents the convergence investigation of the nondimensional natural frequency {\mathrm{\Omega }}_1 of the SSSS homogenous plate. A classical plate (\mu =0\mathrm{n}\mathrm{m}) and nanoplate (\mu =2\mathrm{n}\mathrm{m}) were exposed to a sinusoidally distributed load with an average thickness of a/h=10. To evaluate the convergence rate of the proposed IGA approach, seven alternative finite element meshes were examined for each plate example. As listed in Tab.1, although quick convergence of the analysis is documented for polynomial orders p=3 and p=4, solutions using p=2 exhibit a comparatively moderate convergence rate toward analytical solutions as described by Ref. [42]. Based on the aforementioned findings, an 11 × 11 cubic (p=4) NURBS element mesh was adequate for all the analysis instances. This mesh is used throughout this study, unless otherwise specified.

Next, Tab.2 provides the dimensionless deflection w_1 and dimensionless natural frequency {\mathrm{\Omega }}_1 of a square functionally graded (FG) nanoplate with five different values of the power-law index n_{\mathit{\text{z}}}=0,0.5,2.5,5.5, and 10.5. In the case of plates (\mu =0\mathrm{n} \mathrm{m}) and nanoplates (\mu =2\mathrm{n} \mathrm{m}), it is clear that the outcomes of the current technique are quite similar to those produced by Navier’s solution [42], with only very slight inaccuracies.

In subsequent comparisons, the accuracies of the dimensionless natural frequency {\mathrm{\Omega }}_1 of the isotropic skew plates without an elastic foundation for SSSS and CCCC boundary conditions are listed in Tab.2. These findings can be realized using the model of the study if parameters n_{\mathit{\text{z}}}=0,\lambda =0 ,\alpha =0,\mu =0 , and \beta =0 are provided. They were published by Liew et al. [49], utilizing the Mindlin theory and Rayleigh–Ritz technique. The findings obtained using the model presented in this study correspond with those published by Liew et al. [49], with the highest error corresponding to 2.4949%. When angle \phi is large, the natural frequency in this study is higher when the ratio of a/h=10. Owing to the significant advancements made in the shear strain theory in this study, the deformation energy of the plate reached its fullest expression.

The results in Tab.3 are compared with those of Ebrahimi and Habibi [40] using the finite element method combined with TSDPM. The results for validating dimensionless natural frequency {\hat{\mathrm{\Omega }}}_1 for a rectangular plate composed of a porous material with an asymmetric distribution are listed in Tab.4. This table lists the findings of the investigation of the conditions of liquid saturation. A square plate composed of porous material with an unequal and asymmetrical distribution (E=69\mathrm{G}\mathrm{P}\mathrm{a}; \rho =2260\mathrm{k}\mathrm{g}/{\mathrm{m}}^3;v= 0,25) links SSSS and CCCC with plate size ratios of a /h =5,20. In conditions when the plate size ratio a/ h and Skempton coefficients diverge, it is evident that the findings of this study are congruent with those of Ebrahim and Habibi (the largest error between the two methods is only 1.1872%).

The calculation program and IGA method incorporated into the study are trustworthy as verified by the findings presented above with respect to the issues of static bending and free vibration of the skew plate.

4.2 Static problem

The impact of the variables on the dimensionless deflection w_1 at the mid-plate and dimensionless natural frequency {\text{Ω}}_1 of FGFP skew nanoplates are numerically described in tables and graphs in this section.

The following values were predetermined for each of the survey’s input parameters: a=10\mathrm{n}\mathrm{m},b/a=1, h=1\mathrm{n} \mathrm{m},n_{\mathit{\text{z}}}= 0.5,\lambda =0.2,\alpha =0.3 ,K_{\mathrm{w}}= 20,K_{\mathrm{s}}= 2,\mu =1 \mathrm{n}\mathrm{m}, \varphi = 15°, β = 0.3; unless otherwise specified. The load applied to the plate was uniformly distributed and had a magnitude of q_0.

First, the dimensionless deflection w_1 and dimensionless natural frequency {\text{Ω}}_1 of SSSS FGFP skew nanoplates rest on elastic foundations with three different types of porosities when the Skempton coefficient index is changed with five equally spaced change values of \beta = 0,0.2,0.4,0.6, and 0.8; and the nonlocal factor varies from \mu =0,1,2\mathrm{n} \mathrm{m}. The angle of the plate’s deviation is calculated using values of \varphi = 0°, 15°, 30°, 45° listed in Tab.4 and Tab.5. It is clear from Tab.5 that an increase in the Skempton coefficient index results in a slight reduction in the dimensionless deflection w_1. A potential reason for this is that an increase in coefficient \beta directly affects the material stiffness matrix, which results in an increase in the overall stiffness of the plate. In contrast to the increase in the coefficient \beta, the increase in the nonlocal coefficient \mu decreases the overall stiffness of the structure, thereby lowering the dimensionless deflection w_1. For the skew plate, as the deviation angle \phi increases, the value of dimensionless deflection w_1 decreases, indicating that the stiffness of the plate increases as it grows. This has important implications in the process of calculating and designing FGFP skew nanoplates. For the first dimensionless natural frequency {\text{Ω}}_1 listed in Tab.6, the findings are identical to those for the value of dimensionless deflection w_1, but the increasing or decreasing direction is reversed.

Following this, the displacement value of the dimensionless deflection w_1 and value of the dimensionless natural frequency {\text{Ω}}_1 due to the impact of the porosity coefficient and various porosity types of the skew-nanoplate are presented in Tab.7 and Tab.8. It is evident that an increase in the porosity coefficient results in a significant increase in the dimensionless deflection w_1 as well as a modest increase in the dimensionless natural frequency {\text{Ω}}_1. Given that an increase in the coefficient of porosity results in a reduction in elastic modulus E and mass density \rho, this phenomenon can be explained by the fact that the displacement value w_1 experiences a significant increase when the coefficient of porosity is high. Despite the fact that the frequency value {\text{Ω}}_1 is dependent on E and \rho, in these circumstances, the value of the dimensionless natural frequency {\text{Ω}}_1 increases by a small amount. FGFP 1 is a type of porosity, which was investigated. It produced the least displacement, whereas FGFP 2, which was also investigated, produced the greatest displacement. By gaining an understanding of this process, engineers can develop acceptable porosity types and porosity volumes to limit structural vibrations. However, porosity diminishes the structure to a certain extent.

The influence of different boundary conditions was considered for the plotted findings. In this study, the boundary conditions of CCCC, CFCF, CCCS, CSCS, and SSSS were used. In contrast to analytical approaches, such as that of Navier, which are difficult to obtain, numerical methods are used to resolve diverse boundary conditions and asymmetrical structural forms such as skew plates. Other forms of porosity provide findings that are comparable to those shown in the following section, which depicts the plate structure with FGFP I porosity. Fig.3 and Fig.4 illustrate the effect of plate geometry and various boundary conditions on the dimensionless deflection w_1 and natural frequency {\text{Ω}}_1 of the FGFP skew nanoplates resting on elastic foundations. An increase in the length-to-thickness ratio a/ h (Fig.3 and Fig.4) reduces the value of w_1 and raises the value of the dimensionless natural frequency {\mathrm{\Omega }}_1, and this decrease (increase) is significant. At period a/h = 0.5, the values of w_1 and {\mathrm{\Omega }}_1 remain almost unchanged after this survey. This distinction is an advantage of the high-order shear strain theory developed in this study. For various boundary conditions, the growth of the nonlocal coefficient does not always result in an increase in the dimensionless deflection w_1 (Fig.3(b)). For the specific case of this study, only SSSS and CSCS boundaries enhanced the dimensionless deflection w_1 value, whereas CCCC and CCCS boundaries lowered the dimensionless deflection w_1 value to a certain extent. However, when the coefficient is nonlocal, the dimensionless natural frequency {\text{Ω}}_1 (Fig.4(b)) decreases. This is conceivable due to the fact that the nonlocal coefficient influences the overall stiffness matrix and mass matrix as well as the force vector. The findings obtained due to the modification of the nonlocal coefficient provide engineers with crucial data for nanoplate design computation. The increase in the power-law index n_{\mathit{\text{z}}} (Fig.3(c) and 4(c)) raises the value of the dimensionless deflection w_1 and lowers the dimensionless natural frequency {\text{Ω}}_1. As n_{\mathit{\text{z}}} grows, the density of the metal increases, resulting in a decrease in the rigidity of the plate. An increase in the angle of deflection \phi of the skew plate increases the natural frequency and reduces the displacement value w_1 for all boundary conditions (Fig.3(d) and Fig.4(d)).

Next, Fig.5 and Fig.6 clarify the influence of the porosity coefficient \lambda, Skempton coefficient index \beta, and elastic foundation stiffness K_{\mathrm{w}}, K_{\mathrm{s}} on the values of w_1, {\text{Ω}}_1 of nanoplates with diverse boundary conditions. According to the figure, an increase in the elastic foundation stiffness improves the total plate stiffness, resulting in a reduction in displacement w_1 and an increase in the dimensionless natural frequency {\text{Ω}}_1.

Finally, in this section, the first four vibration mode shapes of the FGFP skew nanoplates resting on elastic foundations with FGFP I porosity and SSSS boundary conditions are shown in Fig.7. We observe that when the plate is rectangular, \phi = 0° with an SSSS symmetric boundary. Furthermore, the displacement of the plate shape when oscillating provides a symmetrical shape, where modes 2 and 3 exhibit similar vibration patterns. However, when the plate is deviated at an angle \phi = 30°, even though the SSSS boundary is symmetric, the displacement of the plate shape deviates by the value of the deviation angle \phi and no vibrational modes exhibit a similar shape.

4.3 Dynamic problems

To assess the validity of the dynamic problem, the deflection of the midpoint over time was examined. Consider an FG plate with the following geometric and mechanical characteristics: E_{\mathrm{c}}= 151\mathrm{G} \mathrm{P}\mathrm{a}, ρ_c = 3000 \mathrm{k}\mathrm{g}/{\mathrm{m}}^3,E_{\mathrm{m}}= 70\mathrm{G} \mathrm{P}\mathrm{a}, {\rho }_{\mathrm{m}}= 2700\mathrm{k}\mathrm{g}/{\mathrm{m}}^3, v_{\mathrm{c}}= v_{\mathrm{m}}= 0.3,a=0.2 \mathrm{m},h =0.01 \mathrm{m}. The plate establishes an SSSS boundary condition for a period of time t under a uniformly distributed load q_0= -{10}^6\mathrm{N} /{\mathrm{m}}^2 and then deloads. Using Navier’s solution and third-order shear deformation theory, this problem model can be derived from the model proposed by Reddy [47].

Over time, without dimensions, tf=t\sqrt{E_{\mathrm{m}}/\left({\rho }_{\mathrm{m}}{a}^2\right)} yields a dimensionless deflection W^{\ast \ast }=w{\displaystyle \frac{E_{\mathrm{m}}h}{q_0{a}^2}} as shown in Fig.8. The findings of this study are comparable to those obtained by Reddy [47] for n_{\mathit{\text{z}}}=0,0.2,2; hence, the method proposed in this study for the dynamic response issue is credible.

Consider an FGFP skew-nanoplate subjected to time-varying loads, including step, triangle, sinusoidal, and explosive loads (using the same input data as in Subsection 4.2). These variable loads are theoretically described as follows: \overline{\mathbb{F}}\left(x,y,t\right)= q_{0}P\left(t\right), where F\left(t\right) denotes the typical time-varying loads of step, triangular, sinusoidal, and explosive blast types, as stated in the following formula and shown in Fig.9.

where t_1 signifies the duration of the structure being subjected to force and t denotes the total survey time. A natural frequency of {\omega }_1 is denoted by notation t_1=0.6t; \tilde{\omega }=0.2 {\omega }_1. The coefficient of the explosive load is \gamma =0.2{\omega }_1. The structural drag coefficient \xi was maintained at a constant value of 0.01 during the survey, which lasted for t=5\times {10}^{-9} \mathrm{s}.

In this section, we examine the time-dependent dimensionless deflection response W_1 at the midpoint of the FGFP skew-nanoplate with FGFP 1 porosity and the effects of the following factors: the nonlocal coefficient, porosity coefficient, deviation angle \phi, Skempton coefficient index \beta, stiffness foundation K_{\mathrm{w}}, K_{\mathrm{s}}, and various boundary conditions.

Fig.10 depicts the influence of the nonlocal coefficient \mu \left(\mathrm{n}\mathrm{m}\right) on the dynamic response of the dimensionless deflection W_1 of SSSS-FGFP I skew nanoplate resting on an elastic foundation. During the application of a dynamic external force, the difference in W_1 value between \mu =0\mathrm{n} \mathrm{m} and \mu =2\mathrm{n}\mathrm{m} is approximately 1.6 times as shown in the figure. After force application for t>0.6t_{\mathrm{m}\mathrm{a}\mathrm{x}}, the plate returns to a condition of free oscillation.

Next, the effect of Skempton coefficient index \beta (shown in Fig.11) and porosity coefficient \lambda (shown in Fig.12) in the material’s microstructure on the dynamic response of dimensionless deflection W_1 of the SSSS-FGFP I skew-nanoplates is provided. We can observe that response W_1 will move in opposing directions when coefficients \beta and \lambda are increased. Nevertheless, an increase in the coefficient \beta only leads to a modest drop in W_1, but an increase in the coefficient \lambda causes a considerable increase in W_1. The difference between a porosity plate with \lambda =0.1 and that with \lambda = 0.8 can be up to 1.3 times.

Finally, the influence of the deviation angle \phi and stiffness foundation K_{\mathrm{w}}, K_{\mathrm{s}} on the dynamic response of the dimensionless deflection W_1 of the SSSS-FGFP I skew-nanoplates is shown in Fig.13 and Fig.14. As shown, an increase in either of the coefficients \varphi , K_{\mathrm{w}}, K_{\mathrm{s}} causes stiffening of the plate, which in turn results in a reduction in the plate’s displacement response W_1. By discovering the laws that govern the increase and decrease in W_1, as a result of the influencing elements, structural designers can gain an overview in a timely manner and provide helpful suggestions in designing plans for the future.

4.4 Other plate geometries including circular and elliptical nanoplates

To determine the effectiveness and applicability of the IGA approach, experiments identical to those of the skewed nanoplate are conducted on elliptical nanoplates as shown in Fig.15. The initial data is as follows: R_y/R_x= 2,h=1 \mathrm{n}\mathrm{m}, n_{\mathit{\text{z}}}=1,\lambda =0.2 ,\alpha =0.3,K_{\mathrm{w}}=20,K_{\mathrm{s}}= 2,\mu =1\mathrm{n} \mathrm{m},\beta =0.3. The radius along the x-axis of the ellipse is as follows: R_x= a/ 2=5\mathrm{n}\mathrm{m}. The findings of deflection w_1 at the center of the plate and natural frequency {\mathrm{\Omega }}_1 of the FGFP elliptical nanoplates resting on elastic foundations with respect to various porosities, nonlocal coefficients, and the Skempton coefficient index are shown in Tab.9 and Tab.10, respectively. Similar to the skew-nanoplate, the midplate deflection increases as the nonlocal coefficient increases. This growth is independent of the geometrical characteristics of the plate. For elliptical plates, it can be observed that with a variable radius R_y/R_x when R_x is constant, the center deflection increases as R_y increases. Consequently, the stiffness of the plate decreases as its overall surface area increases. As Skempton coefficient index \beta increases, the natural frequency of the elliptical nanoplate increases.

As a last point of interest, some deflection response findings at the midpoint of the elliptical plate, which is exposed to typical triangle loads, are presented in Fig.16. The observed results are completely consistent with the laws of physics presented for the skew-nanoplate. Additionally, based on the findings of the displacement response and Fig.16(d), it can be observed that when the load is removed, the plate oscillates freely with a decreasing amplitude. This is the result of structural resistance. If all forms of resistance are disregarded, then the plate oscillates indefinitely with a constant amplitude. It is possible to observe the first four vibration mode configurations of an elliptical nanoplate and circular nanoplate by observing Fig.17. A considerable shift in the vibrational mode can be observed when there is a change in the geometry of the plate. Second and third modes of the circular plate are symmetrical and exhibit the same form; however, the mode of the elliptical plate is unique.

5 Conclusions

In this study, the results of static bending, free vibration, and dynamics of FGFP skew-nanoplates and elliptical nanoplates resting on an elastic foundation are detailed. These results were obtained using the four-unknown shear deformation theory and IGA method. Based on the elastic theory, the relationship between stress and strain was computed for Biot porous materials. Design engineers can select and modify the porosity distribution, porosity distribution density, geometrical parameters, elastic foundation parameters, nonlocal coefficient, and material volume control factor of the plate to obtain the desired results in the calculation process. This can aid in the design of skew-nanoplates and elliptical nanoplates and structures using porous materials. The following are some of the most notable findings of this study.

1) When the porosities of the porous plates are saturated with water, the overall stiffness of the structure increases, resulting in a decrease in deflection and an increase in the fundamental natural frequency of the plate.

2) By increasing the plate’s Skempton coefficient \beta and coefficient of elastic foundation stiffness, the overall stiffness of the plate is also significantly enhanced.

In nanostructures, an increase in the nonlocal coefficient \mu leads to a reduction in the stiffness of the structure; therefore, the deflection is increased and natural frequency is decreased. Conversely, when deviation angle \phi increases with the skew-nanoplate, the plate becomes stiffer.

3) An increase in any of the porosity coefficients \lambda and power law index n_{\mathit{\text{z}}} increases deflection. This is due to the fact that when porosities appear, the stiffness of the plate is significantly reduced.

The following are some of the new points obtained from the study:

4) The combination of the theory of four-unknown shear deformation with the IGA method increases data precision.

5) The combination of the theory of nonlocal elasticity of Enrigen and linear poroelasticity theory of Biot aids in expanding and increasing the ability to investigate the behavior of nanostructures.

The bending, free, and forced vibration behaviors of an FGFP nanoplate with an elastic foundation were investigated.

In future studies on static and dynamic analyses of FGFP nanostructures, the numerical findings of this study will serve as a standard for comparison. The computation software developed in this study will be used in the future for research based on this study. The investigation covers a wide range of subjects, including linear and nonlinear static and dynamic stability concerns and nonlinear static and dynamic analysis.

References

[1] Nguyen V P, Anitescu C, Bordas S P A, Rabczuk T. Isogeometric analysis: An overview and computer implementation aspects. Mathematics and Computers in Simulation, 2015, 117: 89–116

[2] Jüttler B, Langer U, Mantzaflaris A, Moore S E, Zulehner W. Geometry + simulation modules: Implementing isogeometric analysis. Proceedings in Applied Mathematics and Mechanics, 2014, 14(1): 961–962

[3] Marussig B, Hughes T J R. A review of trimming in isogeometric analysis: Challenges, data exchange and simulation aspects. Archives of Computational Methods in Engineering, 2018, 25(4): 1059–1127

[4] Hughes T J R, Cottrell J A, Bazilevs Y. Isogeometric analysis: CAD, finite elements, NURBS, exact geometry and mesh refinement. Computer Methods in Applied Mechanics and Engineering, 2005, 194(39−41): 4135–4195

[5] Hughes T J R, Reali A, Sangalli G. Duality and unified analysis of discrete approximations in structural dynamics and wave propagation: Comparison of p-method finite elements with k-method NURBS. Computer Methods in Applied Mechanics and Engineering, 2008, 197(49−50): 4104–4124

[6] CottrellJ AReali ABazilevsYHughesT J R. Isogeometric analysis of structural vibrations. Computer Methods in Applied Mechanics and Engineering, 2006, 195(41−43): 5257−5296

[7] HughesT J RReali ASangalliG. Efficient quadrature for NURBS-based isogeometric analysis. Computer Methods in Applied Mechanics and Engineering, 2010, 199(5−8): 301−313

[8] DörfelM RJüttlerBSimeonB. Adaptive isogeometric analysis by local h-refinement with T-splines. Computer Methods in Applied Mechanics and Engineering, 2010, 199(5−8): 264−275

[9] BuffaASangalli GVazquezR. Isogeometric analysis in electromagnetics: B-splines approximation. Computer Methods in Applied Mechanics and Engineering, 2010, 199(17−20): 1143−1152

[10] Bazilevs Y, Akkerman I. Large eddy simulation of turbulent Taylor−Couette flow using isogeometric analysis and the residual-based variational multiscale method. Journal of Computational Physics, 2010, 229(9): 3402–3414

[11] CohenEMartin TKirbyR MLycheTRiesenfeld R F. Analysis-aware modeling: Understanding quality considerations in modeling for isogeometric analysis. Computer Methods in Applied Mechanics and Engineering, 2010, 199(5−8): 334−356

[12] Valizadeh N, Natarajan S, Gonzalez-Estrada O A, Rabczuk T, Bui T Q, Bordas S P A. NURBS-based finite element analysis of functionally graded plates: Static bending, vibration, buckling and flutter. Composite Structures, 2013, 99: 309–326

[13] Dsouza S M, Varghese T M, Budarapu P R, Natarajan S. A non-intrusive stochastic isogeometric analysis of functionally graded plates with material uncertainty. Axioms, 2020, 9(3): 92

[14] Hu Q, Xia Y, Natarajan S, Zilian A, Hu P, Bordas S P A. Isogeometric analysis of thin Reissner–Mindlin shells: Locking phenomena and B-bar method. Computational Mechanics, 2020, 65(5): 1323–1341

[15] Ha S H, Choi K K, Cho S. Numerical method for shape optimization using T-spline based isogeometric method. Structural and Multidisciplinary Optimization, 2010, 42(3): 417–428

[16] BazilevsYCalo VZhangYHughesT. Isogeometric fluid-structure interaction analysis with applications to arterial blood flow. Computational Mechanics, 2006, 38(4−5): 310−322

[17] Auricchio F, Da Veiga L B, Hughes T J R, Reali A, Sangalli G. Isogeometric collocation methods. Mathematical Models and Methods in Applied Sciences, 2010, 20(11): 2075–2107

[18] Schmidt R, Kiendl J, Bletzinger K U, Wüchner R. Realization of an integrated structural design process: Analysis-suitable geometric modelling and isogeometric analysis. Computing and Visualization in Science, 2010, 13(7): 315–330

[19] BazilevsJ. Isogeometric analysis of turbulence and fluid-structure interaction. Dissertation for the Doctoral Degree. Austin: The University of Texas at Austin, 2006

[20] Wang Y, Wang Z, Xia Z, Hien Poh L. Structural design optimization using isogeometric analysis: A comprehensive review. Computer Modeling in Engineering & Sciences, 2018, 117(3): 455–507

[21] BontinckZCorno JDeGersem HKurzSPelsA SchöpsSWolfFdeFalco C DölzJVázquez RRömerU. Recent advances of isogeometric analysis in computational electromagnetics. 2017, arXiv: 1709.06004

[22] De Lorenzis L, Wriggers P, Hughes T J R. Isogeometric contact: A review. GAMM Mitteilungen, 2014, 37(1): 85–123

[23] Gao J, Xiao M, Zhang Y, Gao L. A comprehensive review of isogeometric topology optimization: Methods, applications and prospects. Chinese Journal of Mechanical Engineering, 2020, 33(1): 87

[24] Ansari R, Norouzzadeh A. Nonlocal and surface effects on the buckling behavior of functionally graded nanoplates: An isogeometric analysis. Physica E: Low-Dimensional System and Nanostructures, 2016, 84: 84–97

[25] Fan F, Lei B, Sahmani S, Safaei B. On the surface elastic-based shear buckling characteristics of functionally graded composite skew nanoplates. Thin-walled Structures, 2020, 154: 106841

[26] Norouzzadeh A, Ansari R. Isogeometric vibration analysis of functionally graded nanoplates with the consideration of nonlocal and surface effects. Thin-walled Structures, 2018, 127: 354–372

[27] Fan F, Xu Y, Sahmani S, Safaei B. Modified couple stress-based geometrically nonlinear oscillations of porous functionally graded microplates using NURBS-based isogeometric approach. Computer Methods in Applied Mechanics and Engineering, 2020, 372: 113400

[28] Qiu J, Sahmani S, Safaei B. On the NURBS-based isogeometric analysis for couple stress-based nonlinear instability of PFGM microplates. Mechanics Based Design of Structures and Machines, 2020, 51(2): 1–25

[29] Pham Q H, Nguyen P C, Tran V K, Nguyen-Thoi T. Isogeometric analysis for free vibration of bidirectional functionally graded plates in the fluid medium. Defence Technology, 2022, 18(8): 1311–1329

[30] Rahmouni F, Elajrami M, Madani K, Campilho R D S G. Isogeometric analysis based on non-uniform rational B-splines technology of stress and failure strength in inter-ply hybrid laminated composite. Journal of Composite Materials, 2022, 56(18): 2921–2932

[31] Luat D T, Van Thom D, Thanh T T, Van Minh P, Van Ke T, Van Vinh P. Mechanical analysis of bi-functionally graded sandwich nanobeams. Advances in Nano Research, 2021, 11(1): 55–71

[32] Yakoubi K, Montassir S, Moustabchir H, Elkhalfi A, Scutaru M L, Vlase S. Vlase S. T-stress evaluation based cracking of pipes using an extended isogeometric analysis (X-IGA). Symmetry, 2022, 14(5): 1065

[33] Peng X, Lian H, Ma Z, Zheng C. Intrinsic extended isogeometric analysis with emphasis on capturing high gradients or singularities. Engineering Analysis with Boundary Elements, 2022, 134: 231–240

[34] Xia Y, Wang H, Zheng G, Shen G, Hu P. Discontinuous Galerkin isogeometric analysis with peridynamic model for crack simulation of shell structure. Computer Methods in Applied Mechanics and Engineering, 2022, 398: 115193

[35] Sun Y, Zhou Y, Ke Z, Tian K, Wang B. Stiffener layout optimization framework by isogeometric analysis-based stiffness spreading method. Computer Methods in Applied Mechanics and Engineering, 2022, 390: 114348

[36] Bombarde D S, Agrawal M, Gautam S S, Nandy A. A locking-free formulation for three-dimensional isogeometric analysis. Materials Today: Proceedings, 2022, 66: 1710–1715

[37] Shi P, Dong C. A refined hyperbolic shear deformation theory for nonlinear bending and vibration isogeometric analysis of laminated composite plates. Thin-walled Structures, 2022, 174: 109031

[38] Wang L, Yuan X, Xiong C, Wu H. A priori error analysis of discontinuous Galerkin isogeometric analysis approximations of Burgers on surface. Computer Methods in Applied Mechanics and Engineering, 2022, 390: 114342

[39] Thai L M, Luat D T, Phung V B, Van Minh P, Van Thom D. Finite element modeling of mechanical behaviors of piezoelectric nanoplates with flexoelectric effects. Archive of Applied Mechanics, 2022, 92(1): 163–182

[40] Ebrahimi F, Habibi S. Deflection and vibration analysis of higher-order shear deformable compositionally graded porous plate. Steel and Composite Structures, 2016, 20(1): 205–225

[41] Tran V K, Pham Q H, Nguyen-Thoi T. A finite element formulation using four-unknown incorporating nonlocal theory for bending and free vibration analysis of functionally graded nanoplates resting on elastic medium foundations. Engineering with Computers, 2022, 38(2): 1465–1490

[42] Sobhy M. A comprehensive study on FGM nanoplates embedded in an elastic medium. Composite Structures, 2015, 134: 966–980

[43] Biot M A. Theory of buckling of a porous slab and its thermoelastic analogy. Journal of Applied Mechanics, 1964, 31(2): 194–198

[44] Pham Q H, Tran V K, Tran T T, Nguyen-Thoi T, Nguyen P C, Pham V D. A nonlocal quasi-3D theory for thermal free vibration analysis of functionally graded material nanoplates resting on elastic foundation. Case Studies in Thermal Engineering, 2021, 26: 101170

[45] Thai Dung N, Minh Thai L, Van Ke T, Thi Huong Huyen T, Van Minh P. Nonlinear static bending analysis of microplates resting on imperfect two-parameter elastic foundations using modified couple stress theory. Comptes Rendus. Mécanique, 2022, 350(G1): 121–141

[46] Anjomshoa A, Shahidi A R, Hassani B, Jomehzadeh E. Finite element buckling analysis of multi-layered graphene sheets on elastic substrate based on nonlocal elasticity theory. Applied Mathematical Modelling, 2014, 38(24): 5934–5955

[47] ReddyJ N. Analysis of functionally graded plates. International Journal for Numerical Methods in Engineering, 2000, 47(1−3): 663−684

[48] Nguyen N T, Hui D, Lee J, Nguyen-Xuan H. An efficient computational approach for size-dependent analysis of functionally graded nanoplates. Computer Methods in Applied Mechanics and Engineering, 2015, 297: 191–218

[49] Liew K M, Xiang Y, Kitipornchai S, Wang C M. Vibration of thick skew plates based on mindlin shear deformation plate theory. Journal of Sound and Vibration, 1993, 168(1): 39–69

RIGHTS & PERMISSIONS

Higher Education Press