Dynamic analysis of non-homogenous varying thickness rectangular plates resting on Pasternak and Winkler foundations

1. Introduction

Reduction in size, weight and varying density of plates had proven to be of good economic advantage in engineering design. These factors result into better strength, bending, vibration, buckling and guide against resonance as compared to uniform thickness homogenous plates. Also, these increases the utilization of non-homogenous varying thickness plates in engineering fields and industries [1]. In the interest of determining the dynamic analysis of non-homogenous varying thickness rectangular plates, [2-6] investigated the study using analytical method of solution.

The study of structures resting on foundations is of great importance to geotechnics and foundation engineers. However, in the model of foundation in engineering, adoption of Winkler foundation is mostly common which suffers challenge of continuous soil pressure [7]. Reason being that, shear deformation normally occurs in Winkler foundation. To put this wrong assumption in order, designers now introduced Pasternak foundation to support the linear elastic Winkler foundation. However, Sakiyama and Huang [8] utilized green function, a numerical approach to investigate analysis of rectangular plate resting on Winkler foundation. In a further research, Amini et al. [9] showed the power of Chebyshev polynomial and Ritz method in evaluating the analysis of functional graded material rectangular plates on elastic foundation. Bahmyari and Rahbar-Ranji [10] applied Galerkin free element method into the vibration analysis of Orthotropic plates of non-uniform thickness resting on Winkler foundation. In another work, Sharma et al. [11] used differential quadrature solution approach to investigate the behavior of rectangular plate on a Winkler foundation.

Research into dynamic behavior of plates dated back as far as 1800, when Euler first postulated mathematical model for membrane. Thereafter, Lagrange formulated plates governing equation [12]. Furthermore, numerous works had been published in the field. Meanwhile, at initial phase of the research, interest was on simple cases where exact solution can be provided. Chakravery [1] pointed-out that, due to lack of fast computing tools, results obtained for simple cases were not precise, complex cases of vibration problems were left out. This may have led to gap in numerical solution for complex cases of vibration problem despite technological advancement and availability of fast computing tools. It is therefore justified to apply various newly discovered method of solutions in obtaining results for complex and simple cases of vibration problem. Based on literatures reviewed, nonlinear governing differential equations cannot be analyzed using analytical method due to complex mathematics involved. Numerical method has an edge tackling the problem, but it is time consuming with huge cost implication. Nevertheless, in the quest to provide symbolic solutions, Ogunjiofor and Nwoji [13] used characteristic polynomial, a semi-analytical method of solutions to analyze the vibration of plates resting on Winkler foundations. In another study, Matsuda and Sakiyama [14] investigated bending analysis of non-uniform thickness rectangular plates resting on elastic foundation. Furthermore, Miandoab et al. [15] used Homotopy analysis method (HAM) in conducting forced vibration analysis on nanoplates. Thereafter, Rahbar and Rostami [16] discussed vibration response of rectangular orthotropic plates using HAM. HAM suffers limitation of auxiliary parameters and initial approximation. Differential transformation method (DTM), a semi analytical method proposed by Zhou [17] for initial ordinary differential equation, which was later extended to partial differential equation [18], is a reliable method of solution, simple to use without discretization, restrictive assumptions, linearization, perturbation and round-off-error. It reduces complex governing equations into an algebraic equation. Therefore, [19,20] applied DTM on initial value problems of partial differential equations.

The research focuses on application of two-dimensional DTM to nonlinear dynamic behavior of non-homogeneous, non-uniform thickness plate resting on elastic Winkler and Pasternak foundations to establish approximate analytical solutions and parametric studies carried out on the controlling parameters.

2. Problem formulation and mathematical analysis

Considering isotropic rectangular plate of varying thickness shown in Figure 1 resting on two-parameter foundations. Clamped edge boundary condition is considered in this study. The length and width are represented as \(a\) and \(b\) respectively. The governing differential equation of a thin varying thickness rectangular plate as given by [21] is;

\[\begin{aligned} \label{equ1} &&\frac{\partial^2}{\partial x^2} \bigg[ -D \bigg(\frac{\partial^2 w }{\partial x^2} + \nu \frac{\partial^2 w}{\partial y^2} \bigg) \bigg] + 2 \frac{\partial^2}{\partial x \partial y} \bigg[ -D (1- \nu ) \bigg( \frac{\partial^2 w}{\partial x \partial y} \bigg) \bigg] + \frac{\partial^2}{\partial y^2} \bigg[ -D \bigg( \frac{\partial^2 w }{\partial y^2} + \nu \frac{\partial^2 w}{\partial x^2} \bigg) \bigg] + k_w w \notag\\ &&- k_s \bigg( \frac{\partial^2 w}{\partial x^2} + \frac{\partial^2 w}{\partial y^2} \bigg) + k_p w^3 = – \rho h \omega^2 + \frac{\rho h^3 \omega^2}{12} \bigg( \frac{\partial^2 w}{\partial x^2} + \frac{\partial^2 w}{\partial y^2} \bigg) \end{aligned} \tag{1}\] where \(D = \frac{Eh^3}{12(1 – \nu ^2)}\) is the flexural rigidity, \(E\) represents the modulus of rigidity, \(\rho\) represents the density, \(h\) is the thickness, \(w(x,y)\) represents the transverse deflection and \(\omega\) represents the natural frequency, Poisson’s ratio \(v\), shear Pasternak \(k_s\), Winkler stiffness \(k_w\) and \(k_p\) represents nonlinear Winkler stiffness coefficient.

Assuming the thickness and density varies linearly in the direction as:

\[\begin{aligned} \label{equ2} h(x,y) = p_o \bigg( 1 + \beta \frac{x}{a} \bigg) \text{ and homogeneity parameter is } \rho = \rho_o \bigg( 1 + \alpha_1 \frac{x}{a} \bigg) \end{aligned} \tag{2}\] where \(\alpha_1\) is the non-homogeneity parameter and, \(\beta\) is the taper parameter. Therefore, flexural rigidity may be written as [22]:

where the flexural rigidity \(D = \frac{Eh^3}{12(1- \nu ^2)}\). Assuming \(\mu = \frac{Eh^3_m}{12(1-v^2)}\) then, \(D = \mu [h(x,y)]^3\). Therefore,

\[\begin{aligned} \label{equ3} D = \mu \bigg[ 1 + \frac{3x \beta}{a} + \frac{3x^2 \beta^2}{a^2} + \frac{x^3 \beta^3}{a^3} \bigg] \end{aligned} \tag{3}\]

Applying Equation (2) and (3) on Equation (1), we have:

\[\begin{aligned} \label{equ4} &&\frac{\partial^2}{\partial x^2} \bigg[ -\mu h^3 \bigg(\frac{\partial^2 w }{\partial x^2} + \nu \frac{\partial^2 w}{\partial y^2} \bigg) \bigg] + 2 \frac{\partial^2}{\partial x \partial y} \bigg[ – \mu h^3 (1 – \nu ) \bigg( \frac{\partial^2 w}{\partial x \partial y} \bigg) \bigg] + \frac{\partial^2}{\partial y^2} \bigg[ -\mu h^3 \bigg( \frac{\partial^2 w }{\partial y^2} + \nu \frac{\partial^2 w}{\partial x^2} \bigg) \bigg] + k_w w \notag\\&& – k_s \bigg( \frac{\partial^2 w}{\partial x^2} + \frac{\partial^2 w}{\partial y^2} \bigg) + k_p w^3 = – \rho h w \omega^2 + \frac{\rho h^3 \omega^2}{12} \bigg( \frac{\partial^2 w}{\partial x^2} + \frac{\partial^2 w}{\partial y^2} \bigg) \end{aligned} \tag{4}\]

Applying product rule on the Equation (4), we obtained: \[\begin{aligned} \label{equ5} && \mu h^3 \bigg( \frac{\partial^4 w}{\partial x^4} + 2\nu \frac{\partial^4 w}{\partial x^2 \partial y^2} + \frac{\partial^4 w}{\partial y^2} \bigg) + 2 \mu \frac{\partial h^3}{\partial x} \bigg( \frac{\partial^3 w}{\partial x^3} + \nu \frac{\partial^3 w}{\partial y^2 \partial x} \bigg) + 2 \mu \frac{\partial h^3}{\partial y} \bigg( \frac{\partial^3 w}{\partial y^3} + \nu \frac{\partial^3 w}{\partial x^2 \partial y} \bigg) + \mu \frac{\partial^2 h^3}{\partial x^2} \bigg( \frac{\partial^2 w}{\partial x^2} \notag\\ \notag && + \nu \frac{\partial^2 w}{\partial y^2} \bigg)+ \mu \frac{\partial^2 h^3}{\partial y^2} \bigg( \frac{\partial^2 w}{\partial y^2} + \nu \frac{\partial^2 w}{\partial x^2} \bigg) + 2 \mu h^3 (1 -\nu ) \frac{\partial^4 w}{\partial x^2 \partial y^2} + 2 \mu (1- \nu ) \frac{\partial h^3}{\partial x} \frac{\partial^3 w}{\partial x \partial y^2} + 2 \mu (1- \nu ) \frac{\partial h^3}{\partial y} \frac{\partial^3 w}{\partial x^2 \partial y} \\ && + 2 \mu (1 – \nu ) \frac{\partial^2 h^3}{\partial x \partial y} \frac{\partial^2 w}{\partial x \partial y}+ k_ww -k_s \bigg( \frac{\partial^2 w}{\partial x^2} + \frac{\partial^2 w}{\partial y^2} \bigg) + k_pw^3 = \rho h w \omega^2 – \frac{\rho h^3 \omega^2}{12} \bigg( \frac{\partial^2 w}{\partial x^2} + \frac{\partial^2 w}{\partial y^2} \bigg). \end{aligned} \tag{5}\]

Applying the non-dimensionless parameters, \(\xi = \dfrac{x}{a}\), \(\eta = \dfrac{y}{b}\), \(k_w = \dfrac{a^4 \bar{k_w}}{\mu}\), \(k_s = \dfrac{G a^2}{\mu}\), \(k_p = \dfrac{a^4 \bar{k_p}}{\mu}\) and \(W = \dfrac{w }{W_{max}}\) on Equation (5), we arrived at:

\[\begin{aligned} \label{equ6} && \mu \bigg[ h_o \bigg( 1 + \beta \frac{ \xi} {a} \bigg) \bigg]^3 \bigg[ b^4 \frac{\partial^4 W}{ \partial \xi^4 } + 2 a^2 b^2 \frac{\partial^4 W }{\partial \xi^2 \partial \eta^2} + a^4 \frac{\partial^4 W }{\partial \eta^4} \bigg] + 2\mu \frac{\partial \bigg( h_o \bigg( 1 + \beta \frac{\xi}{a} \bigg) \bigg)^3 }{\partial \xi } \bigg[ b^4 \frac{\partial^3 W}{\partial \xi^3} + a^2 b^4 \frac{\partial^3 W}{\partial \xi \partial \eta^2} \bigg] \notag\\&& + \mu \frac{\partial \bigg( \bigg( h_o \bigg( 1 + \beta \frac{\xi}{a} \bigg) \bigg)^3\bigg)^2 }{\partial \xi^2 } \bigg[ b^4 \frac{\partial^2 W}{\partial \xi^2} + a^2 b^4 \nu \frac{\partial^2 W}{\partial \xi^2} \bigg] – a^4 b^4 \rho_o h_o \bigg( 1 + \alpha_1 \frac{\xi }{a} \bigg) \bigg( 1 + \beta \frac{\xi }{a} \bigg) \omega^2 [W] + a^4 b^4 k_w W\notag \\ && – k_s \bigg( a^2 b^4 \frac{\partial^2 W }{\partial \xi^2} + a^4 b^2 \frac{\partial^2 W}{\partial \eta^2} \bigg) + a^4 b^4 k_p W^3 + \frac{\rho_o h_o \bigg( 1 + \alpha_1 \frac{\xi }{a} \bigg) \bigg( 1 + \beta \frac{\xi }{a} \bigg) \omega^2}{12} \bigg[ a^2 b^4 \frac{\partial^2 W }{\partial \xi^2} + a^4 b^2 \frac{\partial^2 W }{\partial \eta^2}\bigg] = 0.\notag\\&& \end{aligned} \tag{6}\]

2.1 Boundary conditions

Clamped edge condition is considered in this study which. This is presented as follow [23]:
Clamped edge (C): \[\begin{aligned} \label{equ7} w(x, y) = 0, \ \ \frac{\partial w(x, y)}{\partial x} = 0 \ at \ x = 0 \ or \ x =a \end{aligned} \tag{7}\] \[\begin{aligned} \label{equ8} w(x, y) = 0, \ \ \frac{\partial w(x, y)}{\partial y} = 0 \ at \ y = 0 \ or \ y =a \end{aligned} \tag{8}\]

3. Method of solution: Two-dimensional differential transformation method

Two-dimensional differential transformation method is an improved version of DTM introduced by Zhou [17]. The method can be applied to partial differential equation with ease. The approach is an iterative method. The basic definitions and operational principles are stated as follows: \[\begin{aligned} \label{equ9} W(k,h) = \frac{1}{k! h!} \bigg[ \frac{\partial^{k+ h} w (x, y)}{\partial x^k \partial y^h} \bigg]_{0,0} , T-\text{function} \end{aligned} \tag{9}\] where \(w(x,y)\) is the original function and \(W(k,h)\) is the transformed function.

The differential inverse transform of \(W(k,h)\) is defined as: \[\begin{aligned} \label{equ10} w(x,y) = \sum_{k =0}^{\infty } \sum_{h = 0}^{\infty } W(x, y) x^k y^h. \end{aligned} \tag{10}\]

From Equation (9) and (10) leads to: \[\begin{aligned} \label{equ11} w(x,y) = \sum_{k =0}^{\infty } \sum_{h = 0}^{\infty }\frac{1}{k! h!} \bigg[ \frac{\partial^{k+ h} w (x, y)}{\partial x^k \partial y^h} \bigg]_{0,0} x^k y^h \end{aligned} \tag{11}\]

3.1 Application of DTM to the boundary conditions

Applying the operation properties of two-dimensional DTM mentioned in Table 2 to the condition Equation (7) at \(x= 0\). As expected, for \(n\)-th order differential governing equation, \(n\)-th number of conditions are required for the analysis. Therefore, for Equation (6), two remaining unknown conditions are represented with constant variables which are found later.
Clamped Support condition: \[\begin{aligned} \label{equ12} W(k,h) = 0 \ \Rightarrow \ W(0, h) = 0, \end{aligned} \tag{12}\] \[\begin{aligned} \label{equ13} (k+ 1) W(k+1, h) = 0 \ \Rightarrow \ W(1, h) = 0. \end{aligned} \tag{13}\]

Therefore, to obtain the four conditions required to analyze fourth-order differential equation as stated earlier, two unknowns are introduced more, we have: \[\begin{aligned} \label{equ14} (k+1 )(K + 2) W (K+ 2, h) = a_1 \ \Rightarrow \ W(2, h) =a_1, \end{aligned} \tag{14}\] \[\begin{aligned} \label{equ15} (k+ 1)(k+ 2)(k+ 3) W (k+3, h) =b_1, \ \Rightarrow \ W(3, h) =b_1. \end{aligned} \tag{15}\]

Therefore, the transformed boundary conditions are: \[\begin{aligned} \label{equ16} W(0, h) = 0; \ W(1, h) = 0; \ W(2, h) =a_1; \ W(3, h) = \frac{b_1}{6}. \end{aligned} \tag{16}\]

3.2. Application of two-dimensional differential transformation method to the governing equation under investigation

Applying the operational properties in Table 1, the two-dimensional DTM transform of the governing equation Equation (6) is resolved as:
\(W_{k+ 4, h} = \frac{1}{\mu {{b}^{4}}{{p}_{0}}^{3}\left( k+1 \right)\left( k+2 \right)\left( k+3 \right)\left( k+4 \right)W(k+4,h)} \times \biggr[ \left( \frac{3\mu {{b}^{4}}\beta {{p}_{0}}^{3}\left( k+1 \right){{\left( k+2 \right)}^{2}}\left( k+3 \right)}{a}\right) W(k+3,h) + 2{{a}^{2}}{{b}^{2}}\mu {{p}_{0}}^{3}\left( k+1 \right)\left( k+2 \right)\left( h+1 \right)\left( h+2 \right)W(k+2,h+2)+ 12\mu a{{b}^{2}}\beta {{p}_{0}}^{3}\left( k+1 \right)kW(k+1,2)+ \left( \frac{4\mu {{b}^{2}}{{\beta }^{3}}{{p}_{0}}^{3}\left( k-1 \right)\left( 3{{b}^{2}}+k-2 \right)}{a} \right)W(k-1,2)+ \mu {{a}^{4}}{{p}_{0}}^{3}\left( h+1 \right)\left( h+2 \right)\left( h+3 \right)\left( h+4 \right)W(k,h+4)+ 3\beta \mu {{a}^{3}}{{p}_{0}}^{3}\left( h+1 \right)\left( h+2 \right)\left( h+3 \right)\left( h+4 \right)W(k-1,h+4)+ 3{{\beta }^{2}}\mu {{a}^{2}}{{p}_{0}}^{3}\left( h+1 \right)\left( h+2 \right)\left( h+3 \right)\left( h+4 \right)W(k-2,h+4)+{{\beta }^{3}}\mu a{{p}_{0}}^{3}\left( h+1 \right)\left( h+2 \right)\left( h+3 \right)\left( h+4 \right)W(k-3,h+4)+ 6\mu a{{b}^{4}}\beta {{p}_{0}}^{3}\left( k+1 \right)\left( h+1 \right)\left( h+2 \right)W(k+1,h+2)+\left( 12\mu {{b}^{2}}{{\beta }^{2}}{{p}_{0}}^{3}k\left( 2{{b}^{2}}+k-1 \right) \right)W(k,2)+ \left( -\frac{\left( k+1 \right)\left( k+2 \right)\left( \left( -\frac{1}{12}{{a}^{4}}{{\omega }^{2}}{{\rho }_{0}}-3{{\beta }^{2}}{{k}^{2}}\mu -9{{\beta }^{2}}k\mu -6{{\beta }^{2}}\mu \right){{p}_{0}}^{3}+{{k}_{s}}{{a}^{4}} \right){{b}^{4}}}{{{a}^{2}}} \right)W(k+2,h)+ \left( -\left( h+1 \right)\left( h+2 \right){{b}^{2}}\left( \left( -\frac{1}{12}{{a}^{4}}{{\omega }^{2}}{{\rho }_{0}}-6{{\beta }^{2}}\mu \nu \right){{p}_{0}}^{3}+{{k}_{s}}{{a}^{4}} \right) \right)W(k,h+2)+ {{a}^{4}}{{b}^{4}}{{k}_{p}}\left(\sum\limits_{l=0}^{k}{\sum\limits_{p=0}^{k-l}{\sum\limits_{r=0}^{h}{\sum\limits_{s=0}^{h-r}{W(l,h-r-s)}W(p,r)W(k-l-p,s)}}} \right)- \left( \frac{{{\omega }^{2}}{{\alpha }_{1}}\left( {{a}^{4}}{{p}_{0}}+\frac{1}{12}{{p}_{0}}^{3}{{\beta }^{2}}\left( k-3 \right)\left( k-2 \right) \right){{\rho }_{0}}{{b}^{4}}}{{{a}^{2}}} \right)W(k-2,h)+\left( \frac{1}{4}\frac{{{p}_{0}}^{3}\left( k+1 \right)k\left( 4\mu \left( k+2 \right)\left( k+1 \right){{\beta }^{3}}+{{a}^{4}}\beta {{\omega }^{2}}{{\rho }_{0}}+\frac{1}{3}{{\rho }_{0}}{{\omega }^{2}}{{a}^{4}}{{\alpha }_{1}} \right){{b}^{4}}}{{{a}^{3}}} \right)W(k+1,h)+{{b}^{4}}\Big( -\left( -\frac{1}{2}k\beta {{p}_{0}}^{3}\left( k-1 \right){{\alpha }_{1}}+{{a}^{4}}{{p}_{0}} \right){{\omega }^{2}}{{\rho }_{0}}+{{a}^{4}}{{k}_{w}} \Big)W(k,h)- \frac{{{\omega }^{2}}\left( -\frac{1}{12}{{p}_{0}}^{3}\left( k-1 \right)\left( k-2 \right){{\beta }^{3}}-\frac{1}{4}{{p}_{0}}^{3}{{\alpha }_{1}}\left( k-1 \right)\left( k-2 \right){{\beta }^{2}}+{{p}_{0}}\beta {{a}^{4}}+{{p}_{0}}{{a}^{4}}{{\alpha }_{1}} \right){{\rho }_{0}}{{b}^{4}}}{a}W(k-1,h)+ \left( \frac{1}{12}\frac{\left( \left( h+1 \right){{p}_{0}}^{3}\left( h+2 \right){{\omega }^{2}}{{\rho }_{0}}{{\alpha }_{1}}+3\beta {{a}^{4}}+72{{\beta }^{3}}\mu \nu \right){{b}^{2}}}{a} \right)W(k-1,h+2)+ \frac{1}{4}{{\rho }_{0}}{{p}_{0}}^{3}{{\omega }^{2}}{{a}^{2}}{{b}^{2}}\beta \left( h+2 \right)\left( h+1 \right)\left( \beta +{{\alpha }_{1}} \right)W(k-2,h+2)+\frac{1}{12}{{\rho }_{0}}{{p}_{0}}^{3}{{\omega }^{2}}a{{b}^{2}}{{\beta }^{2}}\left( h+2 \right)\left( h+1 \right)\left( \beta +3{{\alpha }_{1}} \right)W(k-3,h+2)+ \frac{1}{12}{{\rho }_{0}}{{p}_{0}}^{3}{{\omega }^{2}}{{b}^{2}}{{\beta }^{3}}{{\alpha }_{1}}\left( h+1 \right)\left( h+2 \right)W(k-4,h+2) \biggr].\) \[\begin{aligned} \label{equ17} \end{aligned}\]

3.3. The solution

Applying fixed value of \(h\) and varying \(k\) from 0 to 8 on Equation (17). Similarly reverse the process by varying values of \(h\) from 0 to 8 while \(k\) is fixed, the following equations are developed: \[\begin{aligned} {{W}_{4,0}}&=&-5\beta {{a}_{1}}-\frac{1}{40}\beta {{b}_{1}},\label{equ18}\\ {{W}_{5,0}}&=&\frac{6997{{\beta }^{2}}{{a}_{1}}}{10000}+\frac{3{{\beta }^{2}}{{b}_{1}}}{1000}-\frac{{{b}_{1}}}{30}-\frac{{{a}_{1}}}{5}-\frac{{{b}_{1}}{{\omega }^{2}}}{32967033}+\frac{4{{b}_{1}}{{k}_{s}}}{10989011}+\frac{4262{{a}_{1}}{{\omega }^{2}}}{292721}+\frac{8{{a}_{1}}{{k}_{s}}}{10989011}- \frac{91{{a}_{1}}{{k}_{w}}}{2500000},\label{equ19}\\ {{W}_{6,0}}&=&-\frac{{{\beta }^{3}}{{b}_{1}}}{3000}-\frac{2}{15}{{W}_{4,2}}-\frac{4997{{\beta }^{3}}{{a}_{1}}}{60000}-\frac{849\beta {{a}_{1}}{{\omega }^{2}}}{349951}-\frac{103\beta {{a}_{1}}{{k}_{s}}}{13870933}+\frac{80\beta {{a}_{1}}{{k}_{w}}}{10989011}-\frac{{{b}_{1}}{{\omega }^{2}}{{\alpha }_{1}}}{989010989}\notag\\ &&+ \frac{270{{a}_{1}}{{\omega }^{2}}{{\alpha }_{1}}}{556321}+\frac{\beta {{b}_{1}}{{\omega }^{2}}}{164835165}-\frac{6\beta {{b}_{1}}{{k}_{s}}}{54945055}-\frac{49\beta {{b}_{1}}}{300}+\frac{19999999\beta {{a}_{1}}}{999999951},\label{equ20}\\ {{W}_{7,0}}&=&-\frac{431\beta {{a}_{1}}{{\omega }^{2}}{{\alpha }_{1}}}{4145274}+\frac{\beta {{b}_{1}}{{\omega }^{2}}{{\alpha }_{1}}}{11538461539}+\frac{9{{\omega }^{2}}{{a}_{1}}{{k}_{s}}}{594362519}+\frac{{{\omega }^{2}}{{a}_{1}}{{k}_{w}}}{316990701874}+\frac{4{{\beta }^{2}}{{b}_{1}}{{k}_{s}}}{187617261}\notag \\&&- \frac{{{k}_{s}}{{a}_{1}}{{k}_{w}}}{26415891821}-\frac{{{\beta }^{2}}{{b}_{1}}{{\omega }^{2}}}{1357466065}-\frac{{{\omega }^{2}}{{b}_{1}}{{k}_{s}}}{15849535093702}-\frac{13{{b}_{1}}{{k}_{w}}}{15000000}+\frac{{{a}_{1}}{{k}_{s}}^{2}}{1320794591072} \notag \\&&+ \frac{{{b}_{1}}{{k}_{s}}^{2}}{2641589182144}-\frac{{{\omega }^{4}}{{a}_{1}}}{792480056}+\frac{{{\omega }^{4}}{{b}_{1}}}{380388842345313}-\frac{13{{k}_{p}}{{a}_{1}}^{3}}{2500000}-\frac{14\beta {{W}_{4,2}}}{25}+\frac{4685{{\beta }^{4}}{{a}_{1}}}{504888} \notag\\&&+ \frac{28801{{\beta }^{4}}{{b}_{1}}}{806428001}-\frac{2}{21}{{W}_{5,2}}+\frac{35{{\beta }^{2}}{{a}_{1}}{{k}_{s}}}{15161439}-\frac{13{{\beta }^{2}}{{a}_{1}}{{k}_{w}}}{12500000}+\frac{577{{\beta }^{2}}{{a}_{1}}{{\omega }^{2}}}{1850031}+\frac{13{{a}_{1}}{{\omega }^{2}}{{\alpha }_{1}}}{625000}-\frac{{{b}_{1}}}{210} \notag\\&&+ \frac{4{{a}_{1}}{{\omega }^{2}}}{230769231}+\frac{416{{b}_{1}}{{\omega }^{2}}}{1199995}-\frac{4{{b}_{1}}{{k}_{s}}}{230769231}-\frac{13{{a}_{1}}{{k}_{s}}}{62500000}-\frac{554111{{\beta }^{2}}{{a}_{1}}}{242423563}+\frac{2478{{\beta }^{2}}{{b}_{1}}}{117773},\label{equ21}\\ {{W}_{4,1}}&=&-15\beta {{a}_{1}}-\frac{1}{40}\beta {{b}_{1}},\label{equ22}\\ {{W}_{4,2}}&=&-30\beta {{a}_{1}}-\frac{1}{40}\beta {{b}_{1}},\label{equ23}\\ {{W}_{5,2}}&=&\frac{25991{{\beta }^{2}}{{a}_{1}}}{5000}+\frac{3{{\beta }^{2}}{{b}_{1}}}{1000}-\frac{{{b}_{1}}}{5}-3{{a}_{1}}-\frac{{{b}_{1}}{{\omega }^{2}}}{32967033}+\frac{4{{b}_{1}}{{k}_{s}}}{10989011}+\frac{224{{a}_{1}}{{\omega }^{2}}}{15385}+\frac{48{{a}_{1}}{{k}_{s}}}{10989011}- \frac{91{{a}_{1}}{{k}_{w}}}{2500000},\label{equ24}\\ {{W}_{4,3}}&=&-50\beta {{a}_{1}}-\frac{1}{40}\beta {{b}_{1}},\label{equ25}\\ {{W}_{5,3}}&=&\frac{8797{{\beta }^{2}}{{a}_{1}}}{1000}+\frac{3{{\beta }^{2}}{{b}_{1}}}{1000}-\frac{{{b}_{1}}}{3}-7{{a}_{1}}-\frac{{{b}_{1}}{{\omega }^{2}}}{32967033}+\frac{4{{b}_{1}}{{k}_{s}}}{10989011}+\frac{7023{{a}_{1}}{{\omega }^{2}}}{482369}+\frac{80{{a}_{1}}{{k}_{s}}}{10989011}-\frac{91{{a}_{1}}{{k}_{w}}}{2500000},\label{equ26}\\ {{W}_{4,4}}&=&-75\beta {{a}_{1}}-\frac{1}{40}\beta {{b}_{1}},\label{equ27}\\ {{W}_{5,4}}&=&\frac{26591{{\beta }^{2}}{{a}_{1}}}{2000}+\frac{3{{\beta }^{2}}{{b}_{1}}}{1000}-\frac{{{b}_{1}}}{2}-14{{a}_{1}}-\frac{{{b}_{1}}{{\omega }^{2}}}{32967033}+\frac{4{{b}_{1}}{{k}_{s}}}{10989011}+\frac{1390{{a}_{1}}{{\omega }^{2}}}{95473}+\frac{131{{a}_{1}}{{k}_{s}}}{11996337}- \frac{91{{a}_{1}}{{k}_{w}}}{2500000}.\notag\\&&\label{equ28} \end{aligned} \tag{18-28}\]

Applying DTM principle, we have the solutions of Equation (6) as: \[\begin{aligned} \label{equ29} W(x,y)=\sum\limits_{j=0}^{7}{\sum\limits_{l=0}^{6}{w(j,l){{x}^{j}}{{y}^{l}}}}. \end{aligned} \tag{29}\]

To validate the analytical solutions obtained using two-dimensional differential transformation method. The controlling parameters in the model are assigned value of zeros. Therefore, substitute \({{\alpha }_{1}}\text{ = 0,}\ {{k}_{w}}\text{= 0, }{{k}_{p}}\text{= 0,}\ {{k}_{s}}\text{=0, }\beta \text{ = 0}\) into Equation (29), we obtained:
\(w(x,y)=\frac{{{b}_{1}}{{x}^{3}}{{y}^{6}}}{6}+\frac{{{b}_{1}}{{x}^{3}}{{y}^{3}}}{6}+\frac{{{b}_{1}}{{x}^{3}}y}{6}+{{a}_{1}}x{{y}^{3}}+{{a}_{1}}xy+\frac{{{b}_{1}}{{x}^{3}}{{y}^{4}}}{6}+{{a}_{1}}x{{y}^{5}}+{{a}_{1}}x{{y}^{2}}+\frac{{{b}_{1}}{{x}^{3}}{{y}^{2}}}{6}+\frac{{{b}_{1}}{{x}^{3}}{{y}^{5}}}{6}+ {{a}_{1}}x+\frac{{{b}_{1}}{{x}^{3}}}{6}{{a}_{1}}x{{y}^{6}}+{{a}_{1}}x{{y}^{4}}+\left( -\frac{{{b}_{1}}}{30}-\frac{{{a}_{1}}}{5}-\frac{{{b}_{1}}{{\omega }^{2}}}{32967033}+\frac{4262{{a}_{1}}{{\omega }^{2}}}{292721} \right){{x}^{5}}+ \left( -\frac{{{b}_{1}}}{10}-{{a}_{1}}-\frac{{{b}_{1}}{{\omega }^{2}}}{32967033}+\frac{5824{{a}_{1}}{{\omega }^{2}}}{400005} \right){{x}^{5}}y+\left( \frac{-{{b}_{1}}}{5}-3{{a}_{1}}-\frac{{{b}_{1}}{{\omega }^{2}}}{32967033}+\frac{224{{a}_{1}}{{\omega }^{2}}}{15385} \right){{x}^{5}}{{y}^{2}}+ \left( \frac{-{{b}_{1}}}{3}-7{{a}_{1}}-\frac{{{b}_{1}}{{\omega }^{2}}}{32967033}+\frac{7023{{a}_{1}}{{\omega }^{2}}}{482369} \right){{x}^{5}}{{y}^{3}}+\left( \frac{-{{b}_{1}}}{2}-14{{a}_{1}}-\frac{{{b}_{1}}{{\omega }^{2}}}{32967033}+\frac{1390{{a}_{1}}{{\omega }^{2}}}{95473} \right){{x}^{5}}{{y}^{4}}+ \left( -\frac{7{{b}_{1}}}{10}-\frac{126{{a}_{1}}}{5}-\frac{{{b}_{1}}{{\omega }^{2}}}{32967033}+\frac{2920{{a}_{1}}{{\omega }^{2}}}{200567} \right){{x}^{5}}{{y}^{5}}+\left( -\frac{14{{b}_{1}}}{15}-42{{a}_{1}}-\frac{{{b}_{1}}{{\omega }^{2}}}{32967033}+\frac{1603{{a}_{1}}{{\omega }^{2}}}{110109} \right){{x}^{5}}{{y}^{6}}+ \left( -\frac{1078{{a}_{1}}{{\omega }^{2}}}{777433}+\frac{416{{b}_{1}}{{\omega }^{2}}}{1199985}+\frac{{{\omega }^{4}}{{b}_{1}}}{380388842345313}-\frac{{{\omega }^{4}}{{a}_{1}}}{792480056}+\frac{{{b}_{1}}}{70}+\frac{2{{a}_{1}}}{7} \right){{x}^{7}} \Big( -\frac{1669{{a}_{1}}{{\omega }^{2}}}{401227}+\frac{485{{b}_{1}}{{\omega }^{2}}}{1398986}+\frac{{{\omega }^{4}}{{b}_{1}}}{380388842345313}-\frac{{{\omega }^{4}}{{a}_{1}}}{792486660}+\frac{{{b}_{1}}}{14}+2{{a}_{1}} \Big){{x}^{7}}y+ \Big( -\frac{2628{{a}_{1}}{{\omega }^{2}}}{315895}+\frac{847{{b}_{1}}{{\omega }^{2}}}{2443086}+\frac{{{\omega }^{4}}{{b}_{1}}}{380388842345313}-\frac{{{\omega }^{4}}{{a}_{1}}}{792496567}+\frac{3{{b}_{1}}}{14}+8{{a}_{1}} \Big){{x}^{7}}{{y}^{2}}+ \left( -\frac{3356{{a}_{1}}{{\omega }^{2}}}{242051}+\frac{412{{b}_{1}}{{\omega }^{2}}}{1188313}+\frac{{{\omega }^{4}}{{b}_{1}}}{380388842345313}-\frac{{{\omega }^{4}}{{a}_{1}}}{792509776}+\frac{{{b}_{1}}}{2}+24{{a}_{1}} \right){{x}^{7}}{{y}^{3}}+ \left( -\frac{3949{{a}_{1}}{{\omega }^{2}}}{189889}+\frac{3042{{b}_{1}}{{\omega }^{2}}}{8773355}+\frac{{{\omega }^{4}}{{b}_{1}}}{380388842345313}-\frac{{{\omega }^{4}}{{a}_{1}}}{792526287}+\frac{769230768{{b}_{1}}}{769230769}+60{{a}_{1}} \right){{x}^{7}}{{y}^{4}}+ \Big( -\frac{10141{{a}_{1}}{{\omega }^{2}}}{348327}+\frac{1045{{b}_{1}}{{\omega }^{2}}}{3013632}+\frac{{{\omega }^{4}}{{b}_{1}}}{380388842345313}-\frac{{{\omega }^{4}}{{a}_{1}}}{792546102}+\frac{179999998{{b}_{1}}}{99999999}+132{{a}_{1}} \Big){{x}^{7}}{{y}^{5}}+ \left( -\frac{7751{{a}_{1}}{{\omega }^{2}}}{199687}+\frac{457{{b}_{1}}{{\omega }^{2}}}{1317808}+\frac{{{\omega }^{4}}{{b}_{1}}}{380388842345313}-\frac{{{\omega }^{4}}{{a}_{1}}}{792569221}+3{{b}_{1}}+264{{a}_{1}} \right){{x}^{7}}{{y}^{6}}.\) \[\begin{aligned} \label{equ30} \end{aligned} \tag{30}\]

Applying the boundary conditions Equation (8) on the series solution. Simultaneous equations are obtained for different step values of \(y\) which may be written in this form: for \(x = 1\) and \(y = \frac{2959}{5000}\), Equation (30) becomes: \[\begin{aligned} \label{equ31} W\left[ x,y \right]=\frac{{{\omega }^{4}}{{a}_{1}}}{22760747346362}+\frac{947{{b}_{1}}{{\omega }^{2}}}{572106}+\frac{1621{{a}_{1}}{{\omega }^{2}}}{279789}+\frac{{{\omega }^{4}}{{b}_{1}}}{79662615694136}+\frac{66131{{b}_{1}}}{52363}+ \frac{184670{{a}_{1}}}{51349}=0, \end{aligned} \tag{31}\] \[\begin{aligned} \label{equ32} \frac{dW\left[ x,y \right]}{dx}=\frac{{{\omega }^{4}}{{a}_{1}}}{3793457891060}+\frac{3242{{b}_{1}}{{\omega }^{2}}}{279789}+\frac{84414{{a}_{1}}{{\omega }^{2}}}{2428297}+\frac{{{\omega }^{4}}{{b}_{1}}}{11380373673181}+\frac{180509{{b}_{1}}}{25096}+ \frac{283271{{a}_{1}}}{13723}=0. \end{aligned} \tag{32}\]

Simultaneous Equations (31) and (32) may be written as: \[\left[ \begin{matrix} \left( \dfrac{{{\omega }^{4}}}{22760747346362}+\dfrac{1621{{\omega }^{2}}}{279789}+\dfrac{192992}{53663} \right) & \left( \dfrac{947{{\omega }^{2}}}{572106}+\dfrac{{{\omega }^{4}}}{79662615694136}+\dfrac{8007}{6340} \right) \\ \left( \dfrac{{{\omega }^{4}}}{3793457889621}+\dfrac{84414{{\omega }^{2}}}{2428297}+\dfrac{229581}{11122} \right) & \left( \dfrac{3242{{\omega }^{2}}}{279789}+\dfrac{{{\omega }^{4}}}{11380373671886}+\dfrac{180509}{25096} \right) \\ \end{matrix} \right]\left\{ \begin{matrix} {{a}_{1}} \\ {{b}_{1}} \\ \end{matrix} \right\}=\left\{ \begin{matrix} 0 \\ 0 \\ \end{matrix} \right\} \tag{33}\]

The following characteristic determinant is obtained applying the non-trivial condition: \[\left[ \begin{matrix} \dfrac{{{\omega }^{4}}}{22760747346362}+\dfrac{1621{{\omega }^{2}}}{279789}+\dfrac{192992}{53663} & \dfrac{947{{\omega }^{2}}}{572106}+\dfrac{{{\omega }^{4}}}{79662615694136}+\dfrac{8007}{6340} \\ \dfrac{{{\omega }^{4}}}{3793457889621}+\dfrac{84414{{\omega }^{2}}}{2428297}+\dfrac{229581}{11122} & \dfrac{3242{{\omega }^{2}}}{279789}+\dfrac{{{\omega }^{4}}}{11380373671886}+\dfrac{180509}{25096} \\ \end{matrix} \right] \tag{34}\]

Solving the characteristic determinant results into eigenvalue: \[\begin{aligned} \label{equ35} &&6.59\times {{10}^{-14}}{{\omega }^{4}}+8.69\times {{10}^{-3}}{{\omega }^{2}}+5.39 -1.00\times {{10}^{-22}} \notag\\&& \times \sqrt{( 3.79\times {{10}^{16}}{{\omega }^{8}}+9.99\times {{10}^{28}}{{\omega }^{6}} +6.59\times {{10}^{39}}{{\omega }^{4}} +8.84\times {{10}^{42}}{{\omega }^{2}} +2.93\times {{10}^{45}} )}. \end{aligned} \tag{35}\]

Solving the quadratic Equation (35) results into the natural frequency and the fundamental mode natural frequency is \(\omega =\text{5}\text{.993625349}.\)

3.4. Nonlinear analysis

The Von Kármán principle is adopted for geometric nonlinearity of the equation. The displacement field of the classical thin rectangular plate is modified to account for middle surface displacement components. The strain expressions are expressed as [24]: \[\begin{aligned} {{\varepsilon }_{xx}}&=&\frac{1}{hE}\left( \frac{{{\partial }^{2}}F}{\partial {{y}^{2}}}-\nu \frac{{{\partial }^{2}}F}{\partial {{x}^{2}}} \right),\label{equ36}\\ {{\varepsilon }_{yy}}&=&\frac{1}{hE}\left( \frac{{{\partial }^{2}}F}{\partial {{x}^{2}}}-\nu \frac{{{\partial }^{2}}F}{\partial {{y}^{2}}} \right),\label{equ37}\\ {{\gamma }_{xy}}&=&\frac{-2(1+\nu )}{hE}\frac{{{\partial }^{2}}F}{\partial x\partial y},\label{equ38} \end{aligned} \tag{36–38}\] where \(F\) is the Airy stress function, \(v\) is the Poisson ratio, \(E\) is the Young Modulus and \(h\) is the thickness.

Using these strain expressions (Equation (36) – Equation (38)): \[\begin{aligned} \label{equ39} {{\nabla }^{4}}F=Eh\left( {{\left( \frac{{{\partial }^{2}}W}{\partial x\partial y} \right)}^{2}}-\frac{{{\partial }^{2}}W}{\partial {{x}^{2}}}\frac{{{\partial }^{2}}W}{\partial {{y}^{2}}} \right) \end{aligned} \tag{39}\] \[\begin{aligned} \label{equ40} & {{h}^{3}}\left( \frac{{{\partial }^{4}}W}{\partial {{X}^{4}}}+2{{\lambda }^{2}}\frac{{{\partial }^{4}}W}{\partial {{X}^{2}}\partial {{Y}^{2}}}+{{\lambda }^{4}}\frac{{{\partial }^{4}}W}{\partial {{Y}^{4}}} \right)+2\frac{\partial {{h}^{3}}}{\partial X}\left( \frac{{{\partial }^{3}}W}{\partial {{X}^{3}}}+{{\lambda }^{2}}\frac{{{\partial }^{3}}W}{\partial X\partial {{Y}^{2}}} \right)+\frac{{{\partial }^{2}}{{h}^{3}}}{\partial {{X}^{2}}}\left( \frac{{{\partial }^{2}}W}{\partial {{X}^{2}}}+\nu {{\lambda }^{2}}\frac{{{\partial }^{2}}W}{\partial {{Y}^{2}}} \right)+ \\ \notag & {{k}_{w}}W-{{k}_{s}}\frac{{{\partial }^{2}}W}{\partial {{X}^{2}}}-\frac{{{k}_{s}}}{{{b}^{2}}}\frac{{{\partial }^{2}}W}{\partial {{Y}^{2}}}+{{k}_{p}}{{W}^{3}}=-\rho h\frac{{{\partial }^{2}}W}{\partial {{\tau }^{2}}}+\frac{{{\partial }^{2}}F}{\partial {{Y}^{2}}}\frac{{{\partial }^{2}}W}{\partial {{X}^{2}}}+\frac{{{\partial }^{2}}F}{\partial {{X}^{2}}}\frac{{{\partial }^{2}}W}{\partial {{Y}^{2}}}-2\frac{{{\partial }^{2}}F}{\partial X\partial Y}\frac{{{\partial }^{2}}W}{\partial X\partial Y} \end{aligned} \tag{40}\]

The above Equations (39) and (40) are Von K\(\acute{a}\)rm\(\acute{a}\)n’s plate equations. To obtain an approximate solution for the deflection \(w\), we assume that

• Simply supported: \[\begin{aligned} \label{equ41} w=h\varphi \left( t \right)\cos \left( \frac{\pi }{a} \right)x\cos \left( \frac{\pi }{b} \right)y \end{aligned} \tag{41}\]

• Clamped edge support: \[\begin{aligned} \label{equ42} w=h\varphi \left( t \right){{\cos }^{2}}\left( \frac{\pi }{a} \right)x{{\cos }^{2}}\left( \frac{\pi }{b} \right)y \end{aligned} \tag{42}\]

Applying Equation (41) on Equation (39), we have \[\begin{aligned} \label{equ43} {{\nabla }^{4}}F=\frac{Ehf{{\left( t \right)}^{2}}{{\pi }^{4}}}{{{a}^{2}}{{b}^{2}}}\left( \cos {{\left( \frac{\pi x}{a} \right)}^{2}}\cos {{\left( \frac{\pi y}{b} \right)}^{2}}-\sin {{\left( \frac{\pi x}{a} \right)}^{2}}\sin {{\left( \frac{\pi y}{b} \right)}^{2}} \right) \end{aligned} \tag{43}\] \[\begin{aligned} \label{equ44} {{\nabla }^{4}}F=\frac{1}{2}\frac{Ehf{{\left( t \right)}^{2}}{{\pi }^{4}}}{{{a}^{2}}{{b}^{2}}}\left( \cos \left( \frac{2\pi x}{a} \right)+\cos \left( \frac{2\pi y}{b} \right) \right) \end{aligned} \tag{44}\] \(F\) is the particular integer with complementary function. Particular integer becomes: \[\begin{aligned} \label{equ45} F={{f}_{1}}\left( t \right)\cos \left( \frac{2\pi x}{a} \right)+{{f}_{2}}\left( t \right)\cos \left( \frac{2\pi y}{b} \right) \end{aligned} \tag{45}\] where Equation (44) is: \[\begin{aligned} \label{equ46} \frac{{{\partial }^{4}}F}{\partial {{x}^{4}}}+2\frac{{{\partial }^{4}}F}{\partial {{x}^{2}}\partial {{y}^{2}}}+\frac{{{\partial }^{4}}F}{\partial {{y}^{4}}}=\frac{1}{2}\frac{Ehf{{\left( t \right)}^{2}}{{\pi }^{4}}}{{{a}^{2}}{{b}^{2}}}\left( \cos \left( \frac{2\pi x}{a} \right)+\cos \left( \frac{2\pi y}{b} \right) \right). \end{aligned} \tag{46}\]

Applying Equation (45) on Equation (46) \[\begin{aligned} \label{equ47} \frac{16{{f}_{1}}\left( t \right){{\pi }^{4}}}{{{a}^{4}}}\cos \left( \frac{2\pi x}{a} \right)+\frac{16{{f}_{2}}\left( t \right){{\pi }^{4}}}{{{b}^{4}}}\cos \left( \frac{2\pi y}{b} \right)=\frac{1}{2}\frac{Ehf{{\left( t \right)}^{2}}{{\pi }^{4}}}{{{a}^{2}}{{b}^{2}}}\left( \cos \left( \frac{2\pi x}{a} \right)+\cos \left( \frac{2\pi y}{b} \right) \right) \end{aligned} \tag{47}\]

Equating coefficients of cosine terms in Equation (47), the following are obtained: \[\begin{aligned} \label{equ48} {{f}_{1}}\left( t \right)=\frac{Eh{{\left( f\left( t \right) \right)}^{2}}{{a}^{2}}}{32{{b}^{2}}},{{f}_{2}}\left( t \right)=\frac{Eh{{\left( f\left( t \right) \right)}^{2}}{{b}^{2}}}{32{{a}^{2}}}. \end{aligned} \tag{48}\]

Substituting Equation (48) into Equation (45), we have \[\begin{aligned} \label{equ49} F=\frac{1}{32}Eh{{\left( f\left( t \right) \right)}^{2}}\left( \frac{{{a}^{2}}}{{{b}^{2}}}\cos \left( \frac{2\pi x}{a} \right)+\frac{{{b}^{2}}}{{{a}^{2}}}\cos \left( \frac{2\pi y}{b} \right) \right) \end{aligned} \tag{49}\]

Including the complementary function into Equation (49), one obtains solution for Equation (43) as:

\[\begin{aligned} \label{equ50} F=\frac{1}{2}{{C}_{2}}\left( t \right){{x}^{2}}+\frac{1}{2}{{C}_{1}}\left( t \right){{y}^{2}}+\frac{1}{32}Eh{{\left( f\left( t \right) \right)}^{2}}\left( \frac{{{a}^{2}}}{{{b}^{2}}}\cos \left( \frac{2\pi x}{a} \right)+\frac{{{b}^{2}}}{{{a}^{2}}}\cos \left( \frac{2\pi y}{b} \right) \right). \end{aligned} \tag{50}\]

The in-plane displacement may be written as: \[\ u = \int\limits_{0}^{a}{\left( \frac{1}{hE}\left( \frac{{{\partial }^{2}}F}{\partial {{y}^{2}}} \right)-\nu \frac{{{\partial }^{2}}F}{\partial {{x}^{2}}}-\frac{1}{2}{{\left( \frac{\partial w}{\partial x} \right)}^{2}} \right)}dx,\label{equ51} \tag{51}\] \[ v = \int\limits_{0}^{b}{\left( \frac{1}{hE}\left( \frac{{{\partial }^{2}}F}{\partial {{x}^{2}}} \right)-\nu \frac{{{\partial }^{2}}F}{\partial {{y}^{2}}}-\frac{1}{2}{{\left( \frac{\partial w}{\partial y} \right)}^{2}} \right)}dy.\label{equ52} \tag{52}\]

The boundary condition of the in-plane displacement are: \[\begin{aligned} \label{equ53} u=0, \ and \dfrac{{{\partial }^{2}}F}{\partial x \partial y}=0 \ at \ x=0, a \ and \ 0 \le y \le b \end{aligned} \tag{53}\] \[\begin{aligned} \label{equ54} v=0, \ and \ \dfrac{{{\partial }^{2}}F}{\partial x\partial y}=0 \ at \ y=0, b \ and \ 0 \le x \le a \end{aligned} \tag{54}\] \[\begin{aligned} \label{equ55} {{C}_{1}}\left( t \right)=-\frac{1}{8}\frac{Eh{{\left( f\left( t \right) \right)}^{2}}{{\pi }^{2}}\left( {{a}^{2}}\nu +{{b}^{2}} \right)}{{{a}^{2}}\left( {{\nu }^{2}}-1 \right){{b}^{2}}};\,\,{{C}_{2}}\left( t \right)=-\frac{1}{8}\frac{Eh{{\left( f\left( t \right) \right)}^{2}}{{\pi }^{2}}\left( {{b}^{2}}\nu +{{a}^{2}} \right)}{{{a}^{2}}\left( {{\nu }^{2}}-1 \right){{b}^{2}}} \end{aligned} \tag{55}\]

Hence, the stress function Equation (50) becomes: \[\begin{aligned} \label{equ56} F&=&-\frac{1}{16}\frac{Eh{{\left( f\left( t \right) \right)}^{2}}{{\pi }^{2}}\left( {{b}^{2}}\nu +{{a}^{2}} \right){{x}^{2}}}{{{a}^{2}}\left( {{\nu }^{2}}-1 \right){{b}^{2}}}-\frac{1}{16}\frac{Eh{{\left( f\left( t \right) \right)}^{2}}{{\pi }^{2}}\left( {{a}^{2}}\nu +{{b}^{2}} \right){{y}^{2}}}{{{a}^{2}}\left( {{\nu }^{2}}-1 \right){{b}^{2}}}\notag\\&&+ \frac{1}{32}Eh{{\left( f\left( t \right) \right)}^{2}}\left( \frac{{{a}^{2}}}{{{b}^{2}}}\cos \left( \frac{2\pi x}{a} \right)+\frac{{{b}^{2}}}{{{a}^{2}}}\cos \left( \frac{2\pi y}{b} \right) \right). \end{aligned} \tag{56}\]

Applying Equation (56) and Equation (41) on Equation (40), we have \[\begin{aligned} \label{equ57} \notag L(w,F)&=\left[ 1+\frac{3\beta x}{a}+\frac{3{{\beta }^{2}}{{x}^{2}}}{{{a}^{2}}}+\frac{{{\beta }^{3}}{{x}^{3}}}{{{a}^{3}}} \right] &&\\ \notag & \times \left[ \frac{{{\pi }^{4}}\varphi (t)}{{{a}^{4}}}\cos \left( \frac{\pi x}{a} \right)\cos \left( \frac{\pi y}{b} \right)+\frac{2{{\lambda }^{2}}{{\pi }^{4}}\varphi (t)}{{{a}^{2}}{{b}^{2}}}\cos \left( \frac{\pi x}{a} \right)\cos \left( \frac{\pi y}{b} \right)+ \frac{{{\lambda }^{4}}{{\pi }^{4}}\varphi (t)}{{{b}^{4}}}\cos \left( \frac{\pi x}{a} \right)\cos \left( \frac{\pi y}{b} \right) \right] \\ \notag &+\left[ \frac{6\beta }{a}+\frac{12{{\beta }^{2}}x}{{{a}^{2}}}+\frac{6{{\beta }^{3}}{{x}^{2}}}{{{a}^{3}}} \right]\left( \frac{{{\pi }^{3}}\varphi (t)}{{{a}^{3}}}\sin \left( \frac{\pi x}{a} \right)\cos \left( \frac{\pi y}{b} \right)+\frac{{{\lambda }^{2}}{{\pi }^{3}}\varphi (t)}{a{{b}^{2}}}\sin \left( \frac{\pi x}{a} \right)\cos \left( \frac{\pi y}{b} \right) \right) \\ \notag &+\left[ \frac{6{{\beta }^{2}}}{{{a}^{2}}}+\frac{6{{\beta }^{3}}x}{{{a}^{3}}} \right]\left( -\frac{{{\pi }^{2}}\varphi (t)}{{{a}^{2}}}\cos \left( \frac{\pi x}{a} \right)\cos \left( \frac{\pi y}{b} \right)-\frac{\nu {{\lambda }^{2}}{{\pi }^{2}}\varphi (t)}{{{b}^{2}}}\cos \left( \frac{\pi x}{a} \right)\cos \left( \frac{\pi y}{b} \right) \right) \\ \notag &-\left[ {{\rho }_{o}}{{h}_{o}}+\frac{{{\rho }_{o}}{{h}_{o}}\beta x}{a}+\frac{{{\rho }_{o}}{{h}_{o}}\alpha x}{a}+\frac{{{\rho }_{o}}{{h}_{o}}\alpha {{x}^{2}}\beta }{{{a}^{2}}} \right]\cos \left( \frac{\pi x}{a} \right)\cos \left( \frac{\pi y}{b} \right)\frac{{{\partial }^{2}}W}{\partial {{\tau }^{2}}} \\ \notag &+\frac{{{\pi }^{2}}\varphi (t)}{{{a}^{2}}}\left( -\frac{1}{8}\frac{{{h}^{2}}{{\pi }^{2}}E{{\left( \varphi (t) \right)}^{2}}\left( \nu {{a}^{2}}+{{b}^{2}} \right)}{{{a}^{2}}\left( {{\nu }^{2}}-1 \right){{b}^{2}}}-\frac{1}{8}\frac{{{h}^{2}}{{\pi }^{2}}E{{\left( \varphi (t) \right)}^{2}}}{{{a}^{2}}}\cos \left( \frac{2\pi y}{b} \right) \right)\cos \left( \frac{\pi x}{a} \right)\cos \left( \frac{\pi y}{b} \right) \\ \notag &-\frac{{{\pi }^{2}}\varphi (t)}{{{b}^{2}}}\left( -\frac{1}{8}\frac{{{h}^{2}}{{\pi }^{2}}E{{\left( \varphi (t) \right)}^{2}}\left( {{b}^{2}}\nu +{{a}^{2}} \right)}{{{a}^{2}}\left( {{\nu }^{2}}-1 \right){{b}^{2}}}-\frac{1}{8}\frac{{{h}^{2}}{{\pi }^{2}}E{{\left( \varphi (t) \right)}^{2}}}{{{b}^{2}}}\cos \left( \frac{2\pi x}{a} \right) \right)\cos \left( \frac{\pi x}{a} \right)\cos \left( \frac{\pi y}{b} \right) \\ \notag &+{{k}_{p}}{{\left( \cos \left( \frac{\pi x}{a} \right) \right)}^{3}}{{\left( \cos \left( \frac{\pi y}{b} \right) \right)}^{3}}{{\left( \varphi (t) \right)}^{3}}-{{k}_{w}}\cos \left( \frac{\pi x}{a} \right)\cos \left( \frac{\pi y}{b} \right)\varphi (t) \\ &+{{k}_{s}}\left( -\frac{{{\pi }^{2}}\varphi (t)}{{{a}^{2}}}\cos \left( \frac{\pi x}{a} \right)\cos \left( \frac{\pi y}{b} \right)-\frac{{{\pi }^{2}}\varphi (t)}{{{b}^{2}}}\cos \left( \frac{\pi x}{a} \right)\cos \left( \frac{\pi y}{b} \right) \right)=0. \end{aligned} \tag{57}\]

Meanwhile, \({{h}_{0}}\), \({{\rho }_{o}}\), \(\lambda\) are the thickness of the plate, homogeneity constant at \(x=0({{h}_{n}})\) and aspect ratio respectively.\(\beta\) and \(\alpha\) are the taper constant and the homogeneity constant.

Galerkin equation. \(\int\limits_{0}^{a/2}{\int\limits_{0}^{b/2}{L(w,F)wdxdy=0}}\), we have

\[\begin{aligned} \label{equ58} M{{\ddot{\varphi }}_{s}}(t)+K{{\varphi }_{s}}(t)+V\varphi _{s}^{3}(t)=0 \end{aligned} \tag{58}\] where
\(M =\frac{1}{\biggl( 1536{{\nu }^{2}}-1536 \biggr){{b}^{4}}{{a}^{7}}{{\pi }^{4}}} \biggl[ 12\biggr\{ -8{{a}^{6}}\biggr( -\biggr( \frac{1}{2}b\beta +a \biggr) \biggr( \alpha b+2a \biggr){{\pi }^{2}}+\alpha {{a}^{2}}\beta \biggr)h\biggr( \nu +1 \biggr){{b}^{4}}\rho \biggr( \nu -1 \biggr)\sin \biggr( \frac{b\pi }{a} \biggr) +\biggr( 8{{a}^{7}}h\biggr( \biggr( \alpha +\beta \biggr)a+\alpha b\beta \biggr){{b}^{4}}\rho \biggr( \nu +1 \biggr)\biggr( \nu -1 \biggr)\cos \biggr( \frac{b\pi }{a} \biggr) -8h{{a}^{5}}\biggr( -\frac{1}{2}\biggr( 4{{a}^{2}}+b\biggr( \alpha +\beta \biggr)a+\frac{1}{3}\alpha {{b}^{2}}\beta \biggr)b{{\pi }^{2}}+{{a}^{3}}\biggr( \alpha +\beta \biggr) \biggr)\biggr( \nu +1 \biggr){{b}^{4}}\rho \biggr( \nu -1 \biggr) \biggr)\pi\biggr\} \sin \biggr( \frac{a\pi }{2b} \biggr)b\cos \biggr( \frac{a\pi }{2b} \biggr) +\frac{1}{2}\biggr[ -8{{a}^{6}}\biggr( -\biggr( \frac{1}{2}b\beta +a \biggr)\biggr( \alpha b+2a \biggr){{\pi }^{2}}+\alpha {{a}^{2}}\beta \biggr)h\biggr( \nu +1 \biggr){{b}^{4}}\rho \biggr( \nu -1 \biggr)\sin \biggr( \frac{b\pi }{a} \biggr) +\biggr( 8{{a}^{7}}h\biggr( \biggr( \alpha +\beta \biggr)a+\alpha b\beta \biggr){{b}^{4}}\rho \biggr( \nu +1 \biggr)\biggr( \nu -1 \biggr)\cos \biggr( \frac{b\pi }{a} \biggr) -8h{{a}^{5}}\biggr( -\frac{1}{2}\biggr( 4{{a}^{2}}+b\biggr( \alpha +\beta \biggr)a+\frac{1}{3}\alpha {{b}^{2}}\beta \biggr)b{{\pi}^{2}}+{{a}^{3}}\biggr( \alpha +\beta \biggr) \biggr)\biggr( \nu +1 \biggr){{b}^{4}}\rho \biggr( \nu -1 \biggr) \biggr)\pi \biggr]a\pi \biggr]\).
\(V=\frac{1}{\biggr( 1536{{\nu }^{2}}-1536 \biggr){{b}^{4}}{{a}^{7}}{{\pi }^{4}}}\biggr[ 12\biggr( -2\sin \biggr( \frac{1}{2}\frac{a\pi }{b} \biggr){{a}^{3}}\biggr( \nu +1 \biggr){{b}^{5}}\biggr( a\biggr( E{{\pi }^{4}}{{h}^{2}}-{{k}_{p}}{{a}^{4}}\cos \biggr( \frac{b\pi }{a} \biggr)-4{{k}_{p}}{{a}^{4}} \biggr)\sin \biggr( \frac{b\pi }{a} \biggr) +b\pi \biggr( E{{\pi }^{4}}{{h}^{2}}-3{{k}_{p}}{{a}^{4}} \biggr) \biggr){{\pi }^{2}}\biggr( \nu -1 \biggr){{\biggr( \cos \biggr( \frac{1}{2}\frac{a\pi }{b} \biggr) \biggr)}^{3}} +a\biggr( {{a}^{7}}{{\pi }^{2}}\biggr( \nu -1 \biggr)\biggr( \nu +1 \biggr)\biggr( E{{\pi }^{4}}{{h}^{2}}+3{{b}^{4}}{{k}_{p}} \biggr)\cos \biggr( \frac{b\pi }{a} \biggr) +2{{a}^{3}}\biggr( {{h}^{2}}\biggr( {{a}^{4}}{{\nu }^{2}}-\frac{1}{2}{{b}^{4}}\biggr( {{\nu }^{2}}+1 \biggr) \biggr)E{{\pi }^{4}}+6{{k}_{p}}{{a}^{4}}{{b}^{4}}\biggr( \nu -1 \biggr)\biggr( \nu +1 \biggr) \biggr){{\pi }^{2}} \biggr)\sin \biggr( \frac{b\pi }{a} \biggr) \biggr) \bigg] +\biggr( E{{h}^{2}}\biggr( {{\nu }^{2}}+1 \biggr)\biggr( a-b \biggr)\biggr( a+b \biggr)\biggr( {{a}^{2}}+{{b}^{2}} \biggr){{\pi }^{4}}+9{{k}_{p}}{{a}^{4}}{{b}^{4}}\biggr( \nu -1 \biggr)\biggr( \nu +1 \biggr) \biggr){{a}^{3}}b{{\pi }^{3}}\sin \biggr( \frac{a\pi }{2b} \biggr)b\cos \biggr( \frac{a\pi }{2b} \biggr) +\frac{1}{2}\biggr[ a\biggr( {{a}^{7}}{{\pi }^{2}}\biggr( \nu -1 \biggr)\biggr( \nu +1 \biggr)\biggr( E{{\pi }^{4}}{{h}^{2}}+3{{b}^{4}}{{k}_{p}} \biggr)\cos \biggr( \frac{b\pi }{a} \biggr)+2{{a}^{3}}\biggr( {{h}^{2}}\biggr( {{a}^{4}}{{\nu }^{2}}-\frac{1}{2}{{b}^{4}}\biggr( {{\nu }^{2}}+1 \biggr) \biggr)E{{\pi }^{4}} +6{{k}_{p}}{{a}^{4}}{{b}^{4}}\biggr( \nu -1 \biggr)\biggr( \nu +1 \biggr) \biggr){{\pi }^{2}} \biggr)\sin \biggr( \frac{b\pi }{a} \biggr) +\biggr( E{{h}^{2}}\biggr( {{\nu }^{2}}+1 \biggr)\biggr( a-b \biggr)\biggr( a+b \biggr)\biggr( {{a}^{2}}+{{b}^{2}} \biggr){{\pi }^{4}} +9{{k}_{p}}{{a}^{4}}{{b}^{4}}\biggr( \nu -1 \biggr)\biggr( \nu +1 \biggr) \biggr){{a}^{3}}b{{\pi }^{3}} \biggr]a\pi\).
\(K=\frac{1}{\biggr( 1536{{\nu }^{2}}-1536 \biggr){{b}^{4}}{{a}^{7}}{{\pi }^{4}}}\biggr( 12\biggr( -16a-{{\biggr( 1/2b\beta +a \biggr)}^{3}}{{\biggr( {{a}^{2}}{{\lambda }^{2}}+{{b}^{2}} \biggr)}^{2}}{{\pi }^{4}}+{{a}^{2}}\biggr( \biggr( 3/2{{\lambda }^{4}}{{\beta }^{2}}+{{k}_{s}}{{b}^{2}} \biggr){{a}^{5}}+3/4{{\lambda }^{4}}{{a}^{4}}b{{\beta }^{3}} +\biggr( \biggr( {{k}_{s}}{{b}^{2}}+6\biggr( \nu -1/2 \biggr){{\beta }^{2}}{{\lambda }^{2}} \biggr){{b}^{2}}{{a}^{3}} +3\biggr( \nu -1/2 \biggr){{\beta }^{3}}{{\lambda }^{2}}{{b}^{3}}{{a}^{2}}+3/2a{{b}^{4}}{{\beta }^{2}}+3/4{{b}^{5}}{{\beta }^{3}} \biggr){{\pi }^{2}}+{{k}_{w}}{{a}^{7}}{{b}^{4}} \biggr)\biggr( \nu +1 \biggr)(\nu -1) \biggr) \biggr) {{\pi }^{2}}\sin \biggr( \frac{b\pi }{a} \biggr)+\biggr( 24{{a}^{2}}\biggr( \biggr( {{a}^{2}}{{\lambda }^{2}}-3{{b}^{2}} \biggr){{\biggr( 1/2b\beta +a \biggr)}^{2}}\biggr( {{a}^{2}}{{\lambda }^{2}}+{{b}^{2}} \biggr){{\pi }^{2}}-1/2\biggr( {{\lambda }^{4}}{{a}^{4}}+4\biggr( \nu -1/2 \biggr){{\lambda }^{2}}{{b}^{2}}{{a}^{2}}+{{b}^{4}} \biggr){{a}^{2}}{{\beta }^{2}} \biggr) \beta {{\pi }^{2}}\biggr( \nu +1 \biggr)\biggr( \nu -1 \biggr)\cos \biggr( \frac{b\pi }{a} \biggr)-16\biggr( -\biggr( 1/8{{b}^{2}}{{\beta }^{2}}+1/2ba\beta +{{a}^{2}} \biggr)b\biggr( 1/4b\beta +a \biggr){{\biggr( {{a}^{2}}{{\lambda }^{2}}+{{b}^{2}} \biggr)}^{2}}{{\pi }^{4}} +{{a}^{2}}\biggr( 3/2{{\lambda }^{4}}{{a}^{6}}\beta +{{a}^{5}}{{b}^{3}}{{k}_{s}}-3{{\lambda }^{2}}{{a}^{4}}{{b}^{2}}\beta +{{b}^{3}}\biggr( 6{{\beta }^{2}}{{\lambda }^{2}}\nu +{{k}_{s}}{{b}^{2}} \biggr){{a}^{3}} +3/2{{b}^{4}}\beta \biggr( {{\beta }^{2}}{{\lambda }^{2}}\nu -3 \biggr){{a}^{2}}+6a{{b}^{5}}{{\beta }^{2}}+3/2{{b}^{6}}{{\beta }^{3}} \biggr){{\pi }^{2}} +{{a}^{4}}\biggr( -3/4{{a}^{4}}{{\beta }^{3}}{{\lambda }^{4}}+{{a}^{3}}{{b}^{5}}{{k}_{w}}-3\biggr( \nu -1/2 \biggr){{\beta }^{3}}{{\lambda }^{2}}{{b}^{2}}{{a}^{2}}-3/4{{b}^{4}}{{\beta }^{3}} \biggr)\biggr) \biggr( \nu +1 \biggr)\biggr( \nu -1 \biggr){{\pi }^{2}} \biggr)\pi \sin \biggr( \frac{a\pi }{2b} \biggr)b\cos \biggr( \frac{a\pi }{2b} \biggr) +1/2(-16a(-{{\biggr( 1/2b\beta +a \biggr)}^{3}}{{\biggr( {{a}^{2}}{{\lambda }^{2}}+{{b}^{2}} \biggr)}^{2}}{{\pi }^{4}}+{{a}^{2}}\biggr( 3/2{{\lambda }^{4}}{{\beta }^{2}}+{{k}_{s}}{{b}^{2}} \biggr){{a}^{5}} +3/4{{\lambda }^{4}}{{a}^{4}}b{{\beta }^{3}}+\biggr( {{k}_{s}}{{b}^{2}}+6\biggr( \nu -1/2 \biggr){{\beta }^{2}}{{\lambda }^{2}} \biggr){{b}^{2}}{{a}^{3}}+3\biggr( \nu -1/2 \biggr){{\beta }^{3}}{{\lambda }^{2}}{{b}^{3}}{{a}^{2}} +3/2a{{b}^{4}}{{\beta }^{2}}+3/4{{b}^{5}}{{\beta }^{3}}){{\pi }^{2}}+{{k}_{w}}{{a}^{7}}{{b}^{4}})\biggr( \nu +1 \biggr)\biggr( \nu -1 \biggr){{\pi }^{2}}\sin \biggr( \frac{b\pi }{a} \biggr) +24{{a}^{2}}\biggr( \biggr( {{a}^{2}}{{\lambda }^{2}}-3{{b}^{2}} \biggr){{\biggr( 1/2b\beta +a \biggr)}^{2}}\biggr( {{a}^{2}}{{\lambda }^{2}}+{{b}^{2}} \biggr){{\pi }^{2}} -1/2\biggr( {{\lambda }^{4}}{{a}^{4}}+4\biggr( \nu -1/2 \biggr){{\lambda }^{2}}{{b}^{2}}{{a}^{2}}+{{b}^{4}} \biggr){{a}^{2}}{{\beta }^{2}} \biggr)\beta {{\pi }^{2}}\biggr( \nu +1 \biggr)\biggr( \nu -1 \biggr)\cos \biggr( \frac{b\pi }{a} \biggr) -16\biggr( -\biggr( 1/8{{b}^{2}}{{\beta }^{2}}+1/2ba\beta +{{a}^{2}} \biggr)b\biggr( 1/4b\beta +a \biggr){{\biggr( {{a}^{2}}{{\lambda }^{2}}+{{b}^{2}} \biggr)}^{2}}{{\pi }^{4}} +{{a}^{2}}\biggr( 3/2{{\lambda }^{4}}{{a}^{6}}\beta +{{a}^{5}}{{b}^{3}}{{k}_{s}}-3{{\lambda }^{2}}{{a}^{4}}{{b}^{2}}\beta +{{b}^{3}}\biggr( 6{{\beta }^{2}}{{\lambda }^{2}}\nu +{{k}_{s}}{{b}^{2}} \biggr){{a}^{3}} +3/2{{b}^{4}}\beta \biggr( {{\beta }^{2}}{{\lambda }^{2}}\nu -3 \biggr){{a}^{2}}+6a{{b}^{5}}{{\beta }^{2}}+3/2{{b}^{6}}{{\beta }^{3}} \biggr){{\pi }^{2}} +{{a}^{4}}\biggr( -3/4{{a}^{4}}{{\beta }^{3}}{{\lambda }^{4}}+{{a}^{3}}{{b}^{5}}{{k}_{w}}-3\biggr( \nu -1/2 \biggr){{\beta }^{3}}{{\lambda }^{2}}{{b}^{2}}{{a}^{2}}-3/4{{b}^{4}}{{\beta }^{3}} \biggr) \biggr)\biggr( \nu +1 \biggr)\biggr( \nu -1 \biggr){{\pi }^{2}})\pi )a\pi\).

3.4.1. The initial and boundary conditions

The plate may be subjected to the following boundary conditions:

• Clamped-Clamped support: \[\begin{aligned} \label{equ62} \phi (t)=cosh{{\beta }_{n}}t-cos{{\beta }_{n}}t-\left( \frac{sinh{{\beta }_{n}}L+sin{{\beta }_{n}}L}{cosh{{\beta }_{n}}L-cos{{\beta }_{n}}L} \right)(sinh{{\beta }_{n}}t-sin{{\beta }_{n}}t) \end{aligned} \tag{59}\] where \({{\beta }_{n}}\) are the roots of the equation \(cos{{\beta }_{n}}Lcosh{{\beta }_{n}}L=1\).

  The initial and boundary conditions are \[\begin{aligned} \label{equ63} f(0,t)=a,\,\,\dot{f}(0,t)=0 . f(0,t)={f}'(0,t)=0,\,\,\,\,f(L,t)={f}'(L,t)=0. \end{aligned} \tag{60}\]

  Alternatively, polynomial function of the form Equation (60) can be applied for this type of support system.

  \[\begin{aligned} \label{equ64} \phi (t)=25.20\times ({{t}^{2}}-2{{t}^{3}}+{{t}^{4}}). \end{aligned} \tag{61}\]

• Simply Supported: \[\begin{aligned} \label{equ65} \phi (t)=sin{{\beta }_{n}}t,\,\,\,sin\beta L=0 \text{ implies } {{\beta }_{n}}=\frac{n\pi }{L}. \end{aligned} \tag{62}\]

  The initial and boundary conditions are

  \[\begin{aligned} \label{equ66} f(0,t)=a,\,\,\dot{f}(0,t)=0,\,\,\,\, f(0,t)={f}''(0,t)=0,\,\,\,\,f(L,t)={f}''(L,t)=0. \end{aligned} \tag{63}\]

  Alternatively, polynomial function of the form Equation (63) can be applied for this type of support system. \[\begin{aligned} \label{equ67} \phi (t)=3.20\times (t-2{{t}^{3}}+{{t}^{4}}) \end{aligned} \tag{64}\]

3.5. Determination of natural frequencies

The dynamic response of the structural analysis is carried out under the transformation: \[\begin{aligned} \label{equ68} \tau ={{e}^{iwt}}. \end{aligned} \tag{65}\]

Applying Equation (65) on Equation (58), we have \[\begin{aligned} \label{equ69} M{{\ddot{\varphi }}_{s}}(t)+K{{\varphi }_{s}}(t)-V\varphi _{s}^{3}(t)=0. \end{aligned} \tag{66}\]

In order to find the periodic solution of Equation (66), assume an initial approximation for zero-order deformation as: \[\begin{aligned} \label{equ70} {{\varphi }_{0}}(\tau )=A\cos \tau \end{aligned} \tag{67}\]

Substitute Equation (67) into Equation (66), we have \[\begin{aligned} \label{equ71} -M\omega _{0}^{2}A\cos \tau +KA\cos \tau -V{{A}^{3}}{{\cos }^{3}}\tau =0. \end{aligned} \tag{68}\]

Which gives: \[\begin{aligned} \label{equ72} -M\omega _{0}^{2}A\cos \tau +KA\cos \tau -V{{A}^{3}}\left( \frac{3\cos \tau +\cos 3\tau }{4} \right)=0. \end{aligned} \tag{69}\]

Collect like term, we have \[\begin{aligned} \label{equ73} \left( KA-M\omega _{0}^{2}A-\frac{3V{{A}^{3}}}{4} \right)\cos \tau -\frac{1}{4}V{{A}^{3}}\cos 3\tau =0. \end{aligned} \tag{70}\]

Eliminate secular term, we have \[\begin{aligned} \label{equ74} \left( KA-M\omega _{0}^{2}A-\frac{3V{{A}^{3}}}{4} \right)=0. \end{aligned} \tag{71}\]

Thus, zero-order nonlinear natural frequency becomes: \[\begin{aligned} \label{equ75} {{\omega }_{0}}\approx \sqrt{\frac{K}{M}-\frac{3V{{A}^{2}}}{4M}}. \end{aligned} \tag{72}\]

Therefore, ratio of zero-order nonlinear natural frequency, to the linear frequency: \[\begin{aligned} \label{equ76} \frac{{{\omega }_{0}}}{{{\omega }_{b}}}=\sqrt{1-\frac{3V{{A}^{2}}}{4K}}. \end{aligned} \tag{73}\]

Following the same procedural approach, the first-order nonlinear natural frequency is: \[\begin{aligned} \label{equ77} {{\omega }_{1}}\approx \sqrt{\frac{1}{2}\left\{ \left[ \left( \frac{K}{M} \right)-\left( \frac{3V{{A}^{2}}}{4M} \right) \right]+\sqrt{{{\left[ \left( \frac{K}{M} \right)-\left( \frac{3V{{A}^{2}}}{4M} \right) \right]}^{2}}-\left( \frac{3{{V}^{2}}{{A}^{4}}}{32{{M}^{2}}} \right)} \right\}}. \end{aligned} \tag{74}\]

The ratio of the first-order nonlinear frequency, \(\omega_{1}\) to the linear frequency \(\omega_{b}\) gives: \[\begin{aligned} \label{equ78} \frac{{{\omega }_{1}}}{{{\omega }_{b}}}\approx \sqrt{\frac{1}{2}\left\{ \left[ 1-\left( \frac{3V{{A}^{2}}}{4K} \right) \right]+\sqrt{{{\left[ 1-\left( \frac{3V{{A}^{2}}}{4K} \right) \right]}^{2}}-\left( \frac{3{{V}^{2}}{{A}^{4}}}{32K} \right)} \right\}}. \end{aligned} \tag{75}\]

4. Results and discussion

Numerical solutions obtained using DTM for clamped supported edge and nonlinear analysis are presented here. The material properties for the non-uniform thickness rectangular plates are: Young’s modulus of the plate is \(69.0 GPa\), Poisson’s ratio is \(0.3\), and mass density of the rectangular plate is \(2700 \frac{kg}{m^3}\). The convergence study for the analysis is presented in Table 3. This was carried out by varying the number of iterations for both \(x\) and \(y\) along with each computational time as stated in the table. The study commenced with 5 iterations for \(x\) and 4 for \(y\) gradually till 8, 7. From the results, it is observed that stable convergence is achieved with maximum absolute residual of 8%. Increasing the iterations further increase the precision of the result along with increase in computational time meanwhile at 8,7 no difference in value is obtained. Hence further iterations from this point will increase the computational time and cost with insignificant difference in the established solutions.

The analysis performed is non-dimensionless. Table 4 shows validation of result with cited literature and results are in good agreement.

Influence of the elastic foundation on the fundamental natural frequency of the non-uniform thickness rectangular plates under clamped boundary condition are presented. Tables 5,6,7,8 and 9 present the influence of foundation parameter, taper and homogeneity on the fundamental natural frequency. In this study, consideration is given to:

1. Elastic Winkler type foundation \((k_w = 0,10, 50,100)\)

2. Elastic Pasternak type foundation \((k_s = 0,10,50,100,320)\)

Also, the Figure 2 presents the pictorial view of the relationship between elastic foundation and natural frequency. From the results presented, it is observed that increase in foundation stiffness result in increase in natural frequency of the plate under clamped edge support. Increment is also recorded when plate is resting on two-parameter foundations as presented in Table 5 and Table 6. This is an obvious fact, since the characteristic behavior of the plate is expected to have significant influence on the behavior of the foundation. For the stiffness of both foundation and the plate to be comparable there is need to study the value of stiffness to be adopted for the foundation. Also, the rigidity at the clamped edge support may have contributed to the increase in the natural frequency.

Increase in natural frequency for clamped supported edge condition is not that significant. This is illustrated in Figure 2. This mathematical illustration supports the observation:

\[\begin{aligned} {{\omega }^{2}}=\frac{k}{m}\quad {{f}_{n}}=\frac{1}{2\pi }\sqrt{\frac{k}{m}} \end{aligned} \tag{76}\] where \(k\) is the stiffness, \(m\) is the mass and \(\omega = f_n\) is the natural frequency.

The nonlinear analysis is carried out and the results obtained are illustrated in Figure 3 and 6. The effect of the elastic foundations on variation of the amplitude with ratio of linear to nonlinear natural frequency is studied as shown in Figure 3. From the results, softening nonlinearity is observed. The nonlinear frequency increases as the amplitude of vibration increases. Also, Figure 6 portrays the comparative and combine effect of Winkler and Pasternak foundations on the nonlinear amplitude-frequency response of the rectangular plate. As the combination of Winkler and Pasternak foundations increase, the amplitude of vibration increase. The response is much with inclusion of Pasternak foundation because of the nonlinear foundation parameter.

To consider the influence of the homogeneity on the fundamental natural frequency of the varying thickness rectangular plate, analysis is performed under different values of homogenous parameters adopted which are: \(\alpha = 0\) and \(20\), Taper parameter is assigned zero value to avoid interference of results. From the results obtained and illustrated in Table 6, it is observed that increases homogeneity increase the natural frequency of the rectangular plate under clamped condition. This implies that, as the material properties are being increased, the plates gained more stiffness. The more the plate get stiffened the more it becomes rigid and the natural frequency increases due to increase in stiffness of the plate as a result of increase homogeneity. Tables 5 and 6 also show that, the behavior is general to combine foundations parameter. Figure 4 and 5 illustrate the relationship between homogeneity and natural frequency of combine elastic foundations under clamped support condition.

The observation is justified with the underlined mathematical representation. \[\begin{aligned} natural\ freq.=\frac{\omega {{a}^{2}}}{{{\pi }^{2}}}\sqrt{\frac{\rho h}{D}} \end{aligned} \tag{77}\] where density is \(\rho\), \(D\) is the flexural rigidity, \(h\) is the thickness and \(\Omega\) is the natural frequency.

Tables 8 and 9 illustrate the influence of the taper variation on the natural frequency. Increase in taper ratio along \(x\)-direction of the rectangular plates, means that as the thickness of the plate increase, the plate gain more strength and stiffness. The more the stiffness of the plate, the more that natural frequency of the plate increase. This is an obvious result considering the classical theory of vibration. Figure 7 and 8 show the graphical relationship between the taper increase and the natural frequency of the rectangular plate. This increase is significant based on the results obtained. The negative value of the taper assigned in the Figure 7 and 8 show reduction in the tapering of the plate while the positive values indicate that the thickness are increasing which is shown from zero upward in the graphical illustration respectively. The thread of results is observed for both combine case of foundations and when only elastic Winkler foundation is adopted.

Figure 9 presents the variation on the nonlinear amplitude frequency curve response of the non-uniform rectangular plate. It is observed that as the non-homogeneity values is insignificant with nonlinear natural frequency. The effects of taper parameters variation on the nonlinear amplitude frequency curve response of the rectangular plate are studied and presented in Figure 10. It can be seen that, the nonlinear frequency of the rectangular plate decreases with increase in taper values. Figure 11 illustrates the effect of aspect ratio variation on the nonlinear analysis of the rectangular plate. It is observed that nonlinear frequency of the rectangular plate decrease with the increase aspect ratio values. The aspect ratio has a significant effect on the nonlinear frequency amplitude curve of the plate. It is observed that, an increase in the plate aspect ratio makes the rectangular plate to become stiffer and consequently lowers the nonlinear natural frequency of the plate. Figure 11 shows that values of nonlinear natural frequency have a direct relationship with the aspect ratio values as a result of lowering the softening nonlinearity.

5. Conclusion

The free vibration dynamic behavior of thin rectangular plate resting on elastic Winkler and Pasternak foundation using two-dimensional DTM was investigated. The solutions obtained were validated with literature cited with very precise results. However, the solutions obtained were subsequently used to investigate influence of controlling parameters on elastic foundations. The following observations were recorded.

1. Increase in Winkler foundation leads to increase in natural frequency of the plate. Same increase was recorded in the case of combined Pasternak and Winkler foundations.

2. Increase in taper shows that the natural frequency increases.

3. Increase in non-homogeneity of the plate leads to natural frequency increases.

4. Nonlinear natural frequency ratio increases with increase in elastic foundation parameters.

5. Nonlinear natural frequency decreases with the increase in taper values.

6. It is observed that nonlinear natural frequency decreases with the increase of aspect ratio values.

The present study has revealed the importance of varying thickness, homogeneity on dynamic behaviour of rectangular plates.

Abbreviations: Nomenclature

\(k_w\): Linear Winkler foundation parameters

\(k_s\): Shear modulus of Pasternak foundation

\(\dfrac{d}{dx}\): Differential operator

\(w\): transverse deflection

\(a,b\): Dimension of rectangular plate