Finite element analysis of vibration of an isotropic thin rectangular plate subjected to different boundary conditions under point and uniformly distributed loadings

doi:10.15406/aaoaj.2022.06.00162

eISSN: 2576-4500

Aeronautics and Aerospace Open Access Journal

Research Article Volume 6 Issue 4

Finite element analysis of vibration of an isotropic thin rectangular plate subjected to different boundary conditions under point and uniformly distributed loadings

Sobamowo MG,¹

Verify Captcha

Regret for the inconvenience: we are taking measures to prevent fraudulent form submissions by extractors and page crawlers. Please type the correct Captcha word to see email ID.

Agbelusi CT,¹ Oladosu SA,² Kuku RO²

¹Department of Mechanical Engineering, University of Lagos, Nigeria
²Department of Mechanical Engineering, Lagos State University, Epe Campus, Lagos, Nigeria

Correspondence: Sobamowo MG, Department of Mechanical Engineering, University of Lagos, Nigeria

Received: November 20, 2022 | Published: December 29, 2022

Citation: Sobamowo MG, Agbelusi CT, Oladosu SA, et al. Finite element analysis of vibration of an isotropic thin rectangular plate subjected to different boundary conditions under point and uniformly distributed loadings. Aeron Aero Open Access J. 2022;6(4):196-205. DOI: 10.15406/aaoaj.2022.06.00162

Download PDF

Abstract

In this paper, finite element method is applied to structural analysis of rectangular plate subjected to different boundary conditions. The governing equation of motions for the static analysis of a rectangular plate are developed using Kirchhoff plate bending theory. The verification and the reliability of the finite element methods used is established by comparing the results of an in-house MATLAB code and that of commercial computational solid mechanics code using FLEXPDE as well as exact analytical method. The results of the exact analytical method, MATLAB finite element method and FLEX finite element method show excellent agreements. Also, it was observed that the boundary conditions and the dimensions of the plate have significant influence on the vibration of plates. The results showed that Clamped- Simply Supported - Simply Supported -Simply Supported (CSSS) mixed boundary condition gave the highest deflection for each thickness and loading type while the Clamped-Clamped- Clamped -Simply Supported (CCCS) mixed boundary condition provided the lowest deflection. The FEM algorithm developed in this project shows great agreement with commercial software results. It is inferred from the study that use of the developed finite element codes in MATLAB is an easier, cheaper and quicker approach for the analysis of free vibration of an isotropic rectangular plate.

Nomenclature

a, plate length; b, plate width; E, modulus of elasticity, h, plate thickness; M, Bending moment; P, exciting load per unit area; u in-plane displacements in x-axis; v in-plane displacements in y-axis; v, Poisson’s ratio; w, plate displacement; γ_xy,shear strain; ε_x, in-plane strains in x-axis; ε_y, in-plane strains in y-axis; σ_x, normal stress in x-axis; σ_y, normal stress in y-axis; τ_xy,shear stress; ρ, density

Keywords: structural analysis, rectangular plates, boundary conditions, loading conditions, finite element analysis

Introduction

The studies on static and dynamic behaviours of plates have been presented using various methods.¹ However, vibrational behaviours of complicated shapes and structures with geometry and material nonlinearity as well as nonclassical loading and boundary conditions are very difficult to be accurately described by analytical solutions. Consequently, over the years, the static and dynamic analyses of such structures have been carried out with computational methods. With the aid of meshless Galerkin method, Liew et al.,² numerically analyzed the vibration of corrugated Reissner Mindlin plates while differential quadrature method (DQM) was used by Lu et al.,³ for natural vibration of rectangular plates. Some other researchers have also demonstrated the used of different numerical methods for the free and forced vibrations of plates under differences boundary conditions and loadings. Moon and Choi⁴ adopted Transfer of Dynamic Stiffness Coefficient method for vibration Analysis frame structures. Finite Element-Transfer Stiffness Coefficient Method was utilized by Choi⁵ for the numerical investigations of free vibration of plate structures while Wu and Lu^6,7 demonstrated the use of differential cubature method to investigate the free vibration of arbitrary shaped thick plates. Tanaka et al.,⁸ applied boundary element method to analysis the multi-dimensional structures. Natural vibration of shear deformable general triangular plates was explored by Karunasena and Kitipornchai.⁹ Safaei¹⁰ analyzed the frequency-dependent damped vibrations of multifunctional foam plates. The same author investigated the effect of embedding a porous core on the free vibration behaviour of laminated composite plates.¹¹ Kumar et al.,¹² applied dynamic stiffness method to analyze the free vibration of thin functionally graded rectangular plates. Pratap et al.,¹³ presented vibration analysis of a thin functionally graded plate having an out of plane material inhomogeneity resting on Winkler-Pasternak foundation subjected to different boundary conditions. Ramu and Mohanty¹⁴ studied the natural vibration analysis of rectangular plate structures using finite element method. Liu et al.,¹⁵ presented a unified modeling method for vibration analysis of reinforced functionally graded plate resting on elastic foundation with elastic boundary conditions. Fattahi and Safaei¹⁶ used different nonlocal plate models to study free vibration of graphene sheets in an elastic matrix. Gholami and Ansari¹⁷ explored the nonlinear stability and vibration of pre/post-buckled multilayer functionally graded rectangular plates. In author’s works,^18–20 the same author presented nonlinear vibration of functionally graded graphene platelet-reinforced composite rectangular plates. In this work, finite element method is applied to structural analysis of rectangular plate subjected to different boundary conditions under point and uniformly distribution loads. The governing equation of motions for the static analysis of a rectangular plate are developed using Kirchhoff plate bending theory. The verification and the reliability of the finite element methods used is established by comparing the results of an in-house MATLAB code and that of commercial computational solid mechanics code using FLEXPDE as well as exact analytical method.

Model development

Considering a rectangular plate having thickness, mass density as illustrated in Figure 1. It is assumed that the plate is linearly elastic and the effects of shear deformation as well as that of rotary inertia are neglected. Moreover, the development of the vibration model is based on the Kirchhoff plate bending theory From the basic theory of stress-stress and displacement analysis, the in-plane displacements u (in x-axis) and v (in y-axis) can be expressed as

$u = - z \frac{\partial w}{\partial x}$ $u = - z \frac{\partial w}{\partial x}$ (1)

$v = - z \frac{\partial w}{\partial y}$ $v = - z \frac{\partial w}{\partial y}$ (2)

Figure 1 Free body diagram rectangular plate.

Since the transverse shear deformation is not considered, we can write the in-plane strains as

$ε_{x} = - z \frac{\partial^{2} w}{\partial x^{2}}$ $ε_{x} = - z \frac{\partial^{2} w}{\partial x^{2}}$ (3)

$ε_{y} = - z \frac{\partial^{2} w}{\partial y^{2}}$ $ε_{y} = - z \frac{\partial^{2} w}{\partial y^{2}}$ (4)

$γ_{x y} = - z \frac{\partial^{2} w}{\partial x \partial y}$ $γ_{x y} = - z \frac{\partial^{2} w}{\partial x \partial y}$ (5)

Assuming the plane stress condition for the plate, the constitutive equations for the stresses are given

$σ_{x} = - z D_{x} \frac{\partial^{2} w}{\partial x^{2}}$ $σ_{x} = - z D_{x} \frac{\partial^{2} w}{\partial x^{2}}$ (6)

$σ_{y} = - z D_{y} \frac{\partial^{2} w}{\partial y^{2}}$ $σ_{y} = - z D_{y} \frac{\partial^{2} w}{\partial y^{2}}$ (7)

$τ_{x y} = - z D_{x y} \frac{\partial^{2} w}{\partial x \partial y}$ $τ_{x y} = - z D_{x y} \frac{\partial^{2} w}{\partial x \partial y}$ (8)

Also, the moments as shown in Figure 1 are defined as

$M_{x} = \int_{- h / 2}^{h / 2} σ_{x} z d z$ $M_{x} = \int_{- h / 2}^{h / 2} σ_{x} z d z$ (9)

$M_{y} = \int_{- h / 2}^{h / 2} σ_{y} z d z$ $M_{y} = \int_{- h / 2}^{h / 2} σ_{y} z d z$ (10)

$M_{x y} = \int_{- h / 2}^{h / 2} σ_{x y} z d z$ $M_{x y} = \int_{- h / 2}^{h / 2} σ_{x y} z d z$ (11)

On substitution Eq. (6), (7) and (8) into Eqs. (9), (10) and (11), we have

$M_{x} = - D_{x} \frac{\partial^{2} w}{\partial x^{2}} \int_{- h / 2}^{h / 2} z^{2} d z$ $M_{x} = - D_{x} \frac{\partial^{2} w}{\partial x^{2}} \int_{- h / 2}^{h / 2} z^{2} d z$ (12)

$M_{y} = - D_{y} \frac{\partial^{2} w}{\partial y^{2}} \int_{- h / 2}^{h / 2} z^{2} d z$ $M_{y} = - D_{y} \frac{\partial^{2} w}{\partial y^{2}} \int_{- h / 2}^{h / 2} z^{2} d z$ (13)

$M_{x y} = - D_{x y} \frac{\partial^{2} w}{\partial x \partial y} \int_{- h / 2}^{h / 2} z^{2} d z$ $M_{x y} = - D_{x y} \frac{\partial^{2} w}{\partial x \partial y} \int_{- h / 2}^{h / 2} z^{2} d z$ (14)

Evaluation of Eqs. (12), (13) and (14) gives

$M_{x} = - D_{x} \frac{h^{3}}{12} \frac{\partial^{2} w}{\partial x^{2}}$ $M_{x} = - D_{x} \frac{h^{3}}{12} \frac{\partial^{2} w}{\partial x^{2}}$ (15)

$M_{y} = - D_{y} \frac{h^{3}}{12} \frac{\partial^{2} w}{\partial y^{2}}$ $M_{y} = - D_{y} \frac{h^{3}}{12} \frac{\partial^{2} w}{\partial y^{2}}$ (16)

$M_{x y} = - D_{x y} \frac{h^{3}}{12} \frac{\partial^{2} w}{\partial x \partial y}$ $M_{x y} = - D_{x y} \frac{h^{3}}{12} \frac{\partial^{2} w}{\partial x \partial y}$ (17)

From the free body diagram in Figure 1, we can obtain the equilibrium equation. The moment equilibrium about y- and x-axes and force equilibrium about z-axes yield after neglecting higher order terms

$\frac{\partial M_{x}}{\partial x} + \frac{\partial M_{x y}}{\partial y} - Q_{x} = 0$ $\frac{\partial M_{x}}{\partial x} + \frac{\partial M_{x y}}{\partial y} - Q_{x} = 0$ (18)

$\frac{\partial M_{x y}}{\partial x} + \frac{\partial M_{y}}{\partial y} - Q_{y} = 0$ $\frac{\partial M_{x y}}{\partial x} + \frac{\partial M_{y}}{\partial y} - Q_{y} = 0$ (19)

For a classical plate, the shear forces along x- and y-axes are related to the distributed pressure loading as

$\frac{\partial Q_{x}}{\partial x} + \frac{\partial Q_{y}}{\partial y} + p (x, y) = 0$ $\frac{\partial Q_{x}}{\partial x} + \frac{\partial Q_{y}}{\partial y} + p (x, y) = 0$ (20)

Differentiating Eqs. (18) and (19) with respect to x and y, respectively, we have

$\frac{\partial^{2} M_{x}}{\partial x^{2}} + \frac{\partial^{2} M_{x y}}{\partial x \partial y} = \frac{\partial Q_{x}}{\partial x}$ $\frac{\partial^{2} M_{x}}{\partial x^{2}} + \frac{\partial^{2} M_{x y}}{\partial x \partial y} = \frac{\partial Q_{x}}{\partial x}$ (21)

$\frac{\partial^{2} M_{x y}}{\partial x \partial y} + \frac{\partial^{2} M_{y}}{\partial y^{2}} = \frac{\partial Q_{y}}{\partial y}$ $\frac{\partial^{2} M_{x y}}{\partial x \partial y} + \frac{\partial^{2} M_{y}}{\partial y^{2}} = \frac{\partial Q_{y}}{\partial y}$ (22)

Adding Eqs. (21) and (22), we have

$\frac{\partial^{2} M_{x}}{\partial x^{2}} + 2 \frac{\partial^{2} M_{x y}}{\partial x \partial y} + \frac{\partial^{2} M_{y}}{\partial y^{2}} = \frac{\partial Q_{x}}{\partial x} + \frac{\partial Q_{y}}{\partial y}$ $\frac{\partial^{2} M_{x}}{\partial x^{2}} + 2 \frac{\partial^{2} M_{x y}}{\partial x \partial y} + \frac{\partial^{2} M_{y}}{\partial y^{2}} = \frac{\partial Q_{x}}{\partial x} + \frac{\partial Q_{y}}{\partial y}$ (23)

From Eq. (20)

$\frac{\partial Q_{x}}{\partial x} + \frac{\partial Q_{y}}{\partial y} = - p (x, y)$ $\frac{\partial Q_{x}}{\partial x} + \frac{\partial Q_{y}}{\partial y} = - p (x, y)$ (24)

Substitute Eq. (24) into Eq. (23), we have

$\frac{\partial^{2} M_{x}}{\partial x^{2}} + 2 \frac{\partial^{2} M_{x y}}{\partial x \partial y} + \frac{\partial^{2} M_{y}}{\partial y^{2}} = - p (x, y)$ $\frac{\partial^{2} M_{x}}{\partial x^{2}} + 2 \frac{\partial^{2} M_{x y}}{\partial x \partial y} + \frac{\partial^{2} M_{y}}{\partial y^{2}} = - p (x, y)$ (25)

Eq. (3.25) can be written as

$\frac{\partial^{2} M_{x}}{\partial x^{2}} + 2 \frac{\partial^{2} M_{x y}}{\partial x \partial y} + \frac{\partial^{2} M_{y}}{\partial y^{2}} + p (x, y) = 0$ $\frac{\partial^{2} M_{x}}{\partial x^{2}} + 2 \frac{\partial^{2} M_{x y}}{\partial x \partial y} + \frac{\partial^{2} M_{y}}{\partial y^{2}} + p (x, y) = 0$ (26)

Now, we can substitute Eqs. (3.15), (3.16) and (3.17) into Eq. (3.26),

$\frac{\partial^{2}}{\partial x^{2}} (- D_{x} \frac{h^{3}}{12} \frac{\partial^{2} w}{\partial x^{2}}) + 2 \frac{\partial^{2}}{\partial x \partial y} (- D_{x y} \frac{h^{3}}{12} \frac{\partial^{2} w}{\partial x \partial y}) + \frac{\partial^{2}}{\partial y^{2}} (- D_{y} \frac{h^{3}}{12} \frac{\partial^{2} w}{\partial y^{2}}) + p (x, y) = 0$ $\frac{\partial^{2}}{\partial x^{2}} (- D_{x} \frac{h^{3}}{12} \frac{\partial^{2} w}{\partial x^{2}}) + 2 \frac{\partial^{2}}{\partial x \partial y} (- D_{x y} \frac{h^{3}}{12} \frac{\partial^{2} w}{\partial x \partial y}) + \frac{\partial^{2}}{\partial y^{2}} (- D_{y} \frac{h^{3}}{12} \frac{\partial^{2} w}{\partial y^{2}}) + p (x, y) = 0$ (27)

Evaluation of the above Eq. (27), provides

$\frac{h^{3}}{12} [\frac{\partial^{2}}{\partial x^{2}} (D_{x} \frac{\partial^{2} w}{\partial x^{2}}) + 2 \frac{\partial^{2}}{\partial x \partial y} (D_{x y} \frac{\partial^{2} w}{\partial x \partial y}) + \frac{\partial^{2}}{\partial y^{2}} (D_{y} \frac{\partial^{2} w}{\partial y^{2}})] = p (x, y)$ $\frac{h^{3}}{12} [\frac{\partial^{2}}{\partial x^{2}} (D_{x} \frac{\partial^{2} w}{\partial x^{2}}) + 2 \frac{\partial^{2}}{\partial x \partial y} (D_{x y} \frac{\partial^{2} w}{\partial x \partial y}) + \frac{\partial^{2}}{\partial y^{2}} (D_{y} \frac{\partial^{2} w}{\partial y^{2}})] = p (x, y)$ (28)

For a plate made of homogenous and isotropic material where Dx = Dy = Dz, we have

$\frac{h^{3} D}{12} (\frac{\partial^{4} w}{\partial x^{4}} + 2 \frac{\partial^{4} w}{\partial x^{2} \partial y^{2}} + \frac{\partial^{4} w}{\partial y^{4}}) = p (x, y)$ $\frac{h^{3} D}{12} (\frac{\partial^{4} w}{\partial x^{4}} + 2 \frac{\partial^{4} w}{\partial x^{2} \partial y^{2}} + \frac{\partial^{4} w}{\partial y^{4}}) = p (x, y)$ (29)

where

$D = \frac{E}{1 - ν^{2}}$ $D = \frac{E}{1 - ν^{2}}$

Therefore, Eq. (29) can be written as

$\frac{E h^{3}}{12 (1 - ν^{2})} (\frac{\partial^{4} w}{\partial x^{4}} + 2 \frac{\partial^{4} w}{\partial x^{2} \partial y^{2}} + \frac{\partial^{4} w}{\partial y^{4}}) = p (x, y)$ $\frac{E h^{3}}{12 (1 - ν^{2})} (\frac{\partial^{4} w}{\partial x^{4}} + 2 \frac{\partial^{4} w}{\partial x^{2} \partial y^{2}} + \frac{\partial^{4} w}{\partial y^{4}}) = p (x, y)$ (30)

The expression $\frac{E h^{3}}{12 (1 - ν^{2})}$ $\frac{E h^{3}}{12 (1 - ν^{2})}$ refers to the flexural rigidity of the plate. On denoting $D_{r} = \frac{E h^{3}}{12 (1 - ν^{2})}$ $D_{r} = \frac{E h^{3}}{12 (1 - ν^{2})}$ , we have

$D_{r} (\frac{\partial^{4} w (x, y)}{\partial x^{4}} + 2 \frac{\partial^{4} w (x, y)}{\partial x^{2} \partial y^{2}} + \frac{\partial^{4} w (x, y)}{\partial y^{4}}) = p (x, y)$ $D_{r} (\frac{\partial^{4} w (x, y)}{\partial x^{4}} + 2 \frac{\partial^{4} w (x, y)}{\partial x^{2} \partial y^{2}} + \frac{\partial^{4} w (x, y)}{\partial y^{4}}) = p (x, y)$ (31)

The equation of motion for a plate made of homogenous and isotropic material is

$\frac{\partial^{4} w}{\partial x^{4}} + 2 \frac{\partial^{4} w}{\partial x^{2} \partial y^{2}} + \frac{\partial^{4} w}{\partial y^{4}} = \frac{p (x, y)}{D_{r}}$ $\frac{\partial^{4} w}{\partial x^{4}} + 2 \frac{\partial^{4} w}{\partial x^{2} \partial y^{2}} + \frac{\partial^{4} w}{\partial y^{4}} = \frac{p (x, y)}{D_{r}}$ (32)

For free vibration, we have

$\frac{\partial^{4} w}{\partial x^{4}} + 2 \frac{\partial^{4} w}{\partial x^{2} \partial y^{2}} + \frac{\partial^{4} w}{\partial y^{4}} = 0$ $\frac{\partial^{4} w}{\partial x^{4}} + 2 \frac{\partial^{4} w}{\partial x^{2} \partial y^{2}} + \frac{\partial^{4} w}{\partial y^{4}} = 0$ (33)

Alternatively, we can write

$\nabla^{4} w = 0$ $\nabla^{4} w = 0$ (34)

where the bi-harmonic operator

$\nabla^{4} = \frac{\partial^{4}}{\partial x^{4}} + 2 \frac{\partial^{4}}{\partial x^{2} \partial y^{2}} + \frac{\partial^{4}}{\partial y^{4}}$ $\nabla^{4} = \frac{\partial^{4}}{\partial x^{4}} + 2 \frac{\partial^{4}}{\partial x^{2} \partial y^{2}} + \frac{\partial^{4}}{\partial y^{4}}$ (35)

In the above equations, is the exciting load per unit area, $E, h, ν$ $E, h, ν$ and $ρ (x, y)$ $ρ (x, y)$ are the modulus of elasticity, plate thickness, the Poisson’s ratio, and density, respectively.

The boundary conditions for the static vibration are given

Clamped-Clamped-Clamped-Simply Supported

$\begin{matrix} w (0, y) = \frac{\partial w (0, y)}{\partial x} = 0, w (a, y) = \frac{\partial w (a, y)}{\partial x} = 0 w (x, 0) = \frac{\partial w (x, 0)}{\partial y} = 0, w (x, b) = \frac{\partial^{2} w (x, b)}{\partial y^{2}} = 0 \end{matrix}$ $\begin{array}{l} w (0, y) = \frac{\partial w (0, y)}{\partial x} = 0, w (a, y) = \frac{\partial w (a, y)}{\partial x} = 0 \\ w (x, 0) = \frac{\partial w (x, 0)}{\partial y} = 0, w (x, b) = \frac{\partial^{2} w (x, b)}{\partial y^{2}} = 0 \end{array}$ (36)

Clamped-Simply Supported-Clamped-Simply supported

$\begin{matrix} w (0, y) = \frac{\partial w (0, y)}{\partial x} = 0, w (a, y) = \frac{\partial^{2} w (a, y)}{\partial x^{2}} = 0 w (x, 0) = \frac{\partial w (x, 0)}{\partial y} = 0, w (x, b) = \frac{\partial^{2} w (x, b)}{\partial y^{2}} = 0 \end{matrix}$ $\begin{array}{l} w (0, y) = \frac{\partial w (0, y)}{\partial x} = 0, w (a, y) = \frac{\partial^{2} w (a, y)}{\partial x^{2}} = 0 \\ w (x, 0) = \frac{\partial w (x, 0)}{\partial y} = 0, w (x, b) = \frac{\partial^{2} w (x, b)}{\partial y^{2}} = 0 \end{array}$ (37)

Clamped-Clamped-Simply Supported-Simply Supported

$\begin{matrix} w (0, y) = \frac{\partial w (0, y)}{\partial x} = 0, w (a, y) = \frac{\partial w (a, y)}{\partial x} = 0 w (x, 0) = \frac{\partial^{2} w (x, 0)}{\partial y^{2}} = 0, w (x, b) = \frac{\partial^{2} w (x, b)}{\partial y^{2}} = 0 \end{matrix}$ $\begin{array}{l} w (0, y) = \frac{\partial w (0, y)}{\partial x} = 0, w (a, y) = \frac{\partial w (a, y)}{\partial x} = 0 \\ w (x, 0) = \frac{\partial^{2} w (x, 0)}{\partial y^{2}} = 0, w (x, b) = \frac{\partial^{2} w (x, b)}{\partial y^{2}} = 0 \end{array}$ (38)

Clamped-Simply Supported-Simply-Simply Supported

$\begin{matrix} w (0, y) = \frac{\partial w (0, y)}{\partial x} = 0, w (a, y) = \frac{\partial^{2} w (a, y)}{\partial x^{2}} = 0 w (x, 0) = \frac{\partial^{2} w (x, 0)}{\partial y^{2}} = 0, w (x, b) = \frac{\partial^{2} w (x, b)}{\partial y^{2}} = 0 \end{matrix}$ $\begin{array}{l} w (0, y) = \frac{\partial w (0, y)}{\partial x} = 0, w (a, y) = \frac{\partial^{2} w (a, y)}{\partial x^{2}} = 0 \\ w (x, 0) = \frac{\partial^{2} w (x, 0)}{\partial y^{2}} = 0, w (x, b) = \frac{\partial^{2} w (x, b)}{\partial y^{2}} = 0 \end{array}$ (39)

(Figures 2–5)

Figure 2 Clamped-clamped-clamped-simply supported (CCCS) under loading conditions.

Figure 3 Clamped-simply supported-clamped-simply supported (CSCS) under loading conditions.

Figure 4 Clamped-clamped-simply supported- supported-simply supported (CCSS) under loading conditions.

Figure 5 Clamped-simply supported-simply supported--supported-simply supported (CSSS) under loading conditions.

Finite element analysis to static analysis of rectangular plates

The approximate trial function for the fourth-order partial differential equation is given as

$w (x, y) = 4 \sum v = 1 N_{v} (x, y) W_{v}, v = 1, 2, 3, 4$ $w (x, y) = \sum_{v = 1}^{4} N_{v} (x, y) W_{v}, v = 1, 2, 3, 4$ (40)

where the parametric trial function is taken as

$N (x, y) = α_{1} + α_{2} x + α_{3} y + α_{4} x^{2} + α_{5} x y + α_{6} y^{2},$ $N (x, y) = α_{1} + α_{2} x + α_{3} y + α_{4} x^{2} + α_{5} x y + α_{6} y^{2},$ (41)

The Galerkin weighted residual equation is given by.

$\int Ω N_{v} (x, y) [D_{r} (\frac{\partial^{4} w}{\partial x^{4}} + 2 \frac{\partial^{4} w}{\partial x^{2} \partial y^{2}} + \frac{\partial^{4} w}{\partial y^{4}}) - p] d x d y = 0$ $\int_{Ω} N_{v} (x, y) [D_{r} (\frac{\partial^{4} w}{\partial x^{4}} + 2 \frac{\partial^{4} w}{\partial x^{2} \partial y^{2}} + \frac{\partial^{4} w}{\partial y^{4}}) - p] d x d y = 0$ (42)

After the weak formulation, we have

$\begin{matrix} \int A w D_{r} [\frac{\partial N_{i}^{2}}{\partial x^{2}} \frac{\partial N_{j}^{2}}{\partial x^{2}} + 2 \frac{\partial N_{i}^{2}}{\partial x^{2}} \frac{\partial N_{j}^{2}}{\partial y^{2}} + \frac{\partial N_{i}^{2}}{\partial y^{2}} \frac{\partial N_{j}^{2}}{\partial y^{2}}] - ⎡ ⎢ ⎣ \int s N_{i} θ_{n} d s + \int M_{n} \frac{\partial N_{i}}{\partial n} d s - \int s \frac{\partial M_{n s}}{\partial s} + {[N_{i} M_{n s}]}_{s} - \int A N_{i} p d A ⎤ ⎥ ⎦ = 0, \end{matrix}$ $\begin{array}{l} \int_{A} w D_{r} [\frac{\partial N_{i}^{2}}{\partial x^{2}} \frac{\partial N_{j}^{2}}{\partial x^{2}} + 2 \frac{\partial N_{i}^{2}}{\partial x^{2}} \frac{\partial N_{j}^{2}}{\partial y^{2}} + \frac{\partial N_{i}^{2}}{\partial y^{2}} \frac{\partial N_{j}^{2}}{\partial y^{2}}] \\ - [\int_{s} N_{i} θ_{n} d s + \int M_{n} \frac{\partial N_{i}}{\partial n} d s - \int_{s} \frac{\partial M_{n s}}{\partial s} + {[N_{i} M_{n s}]}_{s} - \int_{A} N_{i} p d A] = 0, \end{array}$ (43)

Which can be written in matrix form as

$[k] [W] = {F}$ $[k] [W] = {F}$ (44)

the local stiffness matrix

$\begin{matrix} [k] = \int A w D_{r} [\frac{\partial N_{i}^{2}}{\partial x^{2}} \frac{\partial N_{j}^{2}}{\partial x^{2}} + 2 \frac{\partial N_{i}^{2}}{\partial x^{2}} \frac{\partial N_{j}^{2}}{\partial y^{2}} + \frac{\partial N_{i}^{2}}{\partial y^{2}} \frac{\partial N_{j}^{2}}{\partial y^{2}}] - ⎡ ⎣ \int s N_{i} θ_{n} d s + \int M_{n} \frac{\partial N_{i}}{\partial n} d s - \int s \frac{\partial M_{n s}}{\partial s} + {[N_{i} M_{n s}]}_{s} ⎤ ⎦ = 0, \end{matrix}$ $\begin{array}{l} [k] = \int_{A} w D_{r} [\frac{\partial N_{i}^{2}}{\partial x^{2}} \frac{\partial N_{j}^{2}}{\partial x^{2}} + 2 \frac{\partial N_{i}^{2}}{\partial x^{2}} \frac{\partial N_{j}^{2}}{\partial y^{2}} + \frac{\partial N_{i}^{2}}{\partial y^{2}} \frac{\partial N_{j}^{2}}{\partial y^{2}}] \\ - [\int_{s} N_{i} θ_{n} d s + \int M_{n} \frac{\partial N_{i}}{\partial n} d s - \int_{s} \frac{\partial M_{n s}}{\partial s} + {[N_{i} M_{n s}]}_{s}] = 0, \end{array}$ (45)

and the local load matrix

${F} = - \int A N_{i} p d A,$ ${F} = - \int_{A} N_{i} p d A,$ (46)

Finite element solutions to the natural frequencies of the rectangular plates under different boundary conditions

From the governing equation in previous section, the finite element solution of the natural frequency of the rectangular plate can be expressed as

$ω^{2} = \frac{\int Ω N_{v} (x, y) [D_{r} (\frac{\partial^{4} w (x, y)}{\partial x^{4}} + 2 \frac{\partial^{4} w (x, y)}{\partial x^{2} \partial y^{2}} + \frac{\partial^{4} w (x, y)}{\partial y^{4}})] d x d y}{\int Ω ρ h N_{v} (x, y) w (x, y) d x d y}$ $ω^{2} = \frac{\int_{Ω} N_{v} (x, y) [D_{r} (\frac{\partial^{4} w (x, y)}{\partial x^{4}} + 2 \frac{\partial^{4} w (x, y)}{\partial x^{2} \partial y^{2}} + \frac{\partial^{4} w (x, y)}{\partial y^{4}})] d x d y}{\int_{Ω} ρ h N_{v} (x, y) w (x, y) d x d y}$ (48)

Which can be well written as

$ω^{2} = \frac{\iint A N_{v} (x, y) [D_{r} (\frac{\partial^{4} w (x, y)}{\partial x^{4}} + 2 \frac{\partial^{4} w (x, y)}{\partial x^{2} \partial y^{2}} + \frac{\partial^{4} w (x, y)}{\partial y^{4}})] d x d y}{\iint A ρ h N_{v} (x, y) w (x, y) d x d y}$ $ω^{2} = \frac{\iint_{A} N_{v} (x, y) [D_{r} (\frac{\partial^{4} w (x, y)}{\partial x^{4}} + 2 \frac{\partial^{4} w (x, y)}{\partial x^{2} \partial y^{2}} + \frac{\partial^{4} w (x, y)}{\partial y^{4}})] d x d y}{\iint_{A} ρ h N_{v} (x, y) w (x, y) d x d y}$ (49)

If we substitute Eq. (42), we have

$ω^{2} = \frac{\iint A N_{v} (x, y) [D_{r} (\frac{\partial^{4}}{\partial x^{4}} (4 \sum v = 1 N_{v} (x, y) W_{v}) + 2 \frac{\partial^{4}}{\partial x^{2} \partial y^{2}} (4 \sum v = 1 N_{v} (x, y) W_{v}) + \frac{\partial^{4}}{\partial y^{4}} (4 \sum v = 1 N_{v} (x, y) W_{v}))] d x d y}{\iint A ρ h N_{v} (x, y) w (x, y) d x d y}$ $ω^{2} = \frac{\iint_{A} N_{v} (x, y) [D_{r} (\frac{\partial^{4}}{\partial x^{4}} (\sum_{v = 1}^{4} N_{v} (x, y) W_{v}) + 2 \frac{\partial^{4}}{\partial x^{2} \partial y^{2}} (\sum_{v = 1}^{4} N_{v} (x, y) W_{v}) + \frac{\partial^{4}}{\partial y^{4}} (\sum_{v = 1}^{4} N_{v} (x, y) W_{v}))] d x d y}{\iint_{A} ρ h N_{v} (x, y) w (x, y) d x d y}$ (50)

Substituting the shape function in Eq. (41) and carrying out the integration, we have

Clamped-clamped-clamped-simply supported

$ω = π^{2} \sqrt{\frac{D_{r}}{ρ h} (\frac{5.1725}{a^{4}} + \frac{2.8050}{a^{2} b^{2}} + \frac{2.4533}{b^{4}})}$ $ω = π^{2} \sqrt{\frac{D_{r}}{ρ h} (\frac{5.1725}{a^{4}} + \frac{2.8050}{a^{2} b^{2}} + \frac{2.4533}{b^{4}})}$ (51)

Clamped-simply supported-clamped-simply supported

$ω = π^{2} \sqrt{\frac{D_{r}}{ρ h} (\frac{1.0016}{a^{4}} + \frac{2.4316}{a^{2} b^{2}} + \frac{5.1747}{b^{4}})}$ $ω = π^{2} \sqrt{\frac{D_{r}}{ρ h} (\frac{1.0016}{a^{4}} + \frac{2.4316}{a^{2} b^{2}} + \frac{5.1747}{b^{4}})}$ (52)

Clamped-clamped-simply supported-simply supported

$ω = π^{2} \sqrt{\frac{D_{r}}{ρ h} (\frac{2.4506}{a^{4}} + \frac{2.6536}{a^{2} b^{2}} + \frac{2.4506}{b^{4}})}$ $ω = π^{2} \sqrt{\frac{D_{r}}{ρ h} (\frac{2.4506}{a^{4}} + \frac{2.6536}{a^{2} b^{2}} + \frac{2.4506}{b^{4}})}$ (53)

Clamped-simply supported-simply-simply supported

$ω = π^{2} \sqrt{\frac{D_{r}}{ρ h} (\frac{5.1747}{a^{4}} + \frac{2.9565}{a^{2} b^{2}} + \frac{5.1747}{b^{4}})}$ $ω = π^{2} \sqrt{\frac{D_{r}}{ρ h} (\frac{5.1747}{a^{4}} + \frac{2.9565}{a^{2} b^{2}} + \frac{5.1747}{b^{4}})}$ (54)

Clamped- clamped-clamped-clamped-supported

Simply supported -simply supported-simply-simply supported

$ω = π^{2} \sqrt{\frac{D_{r}}{ρ h} (\frac{1.0015}{a^{4}} + \frac{2.0008}{a^{2} b^{2}} + \frac{1.0015}{b^{4}})}$ $ω = π^{2} \sqrt{\frac{D_{r}}{ρ h} (\frac{1.0015}{a^{4}} + \frac{2.0008}{a^{2} b^{2}} + \frac{1.0015}{b^{4}})}$ (56)

For the sake of comparison of results, exact analytical solutions for the plate under different boundary conditions are developed. For example, for the Simply Supported -Simply Supported-Simply-Simply Supported plate with point load is given as

$w (x, y) = \infty \sum m = 1 \infty \sum n = 1 \frac{4 p (η, ξ)}{a b} \frac{s i n (\frac{m π}{a}) s i n (\frac{n π}{b})}{ω_{m n}} ⎡ ⎣ \begin{matrix} (s i n (\frac{m π x}{a}) c o s (\frac{m π}{2}) + c o s (\frac{m π x}{a}) s i n (\frac{m π}{2})) \times (s i n (\frac{n π}{b}) c o s (\frac{n π}{2}) + c o s (\frac{n π y}{a}) s i n (\frac{n π}{2})) \end{matrix} ⎤ ⎦$ $w (x, y) = \sum_{m = 1}^{\infty} \sum_{n = 1}^{\infty} \frac{4 p (η, ξ)}{a b} \frac{s i n (\frac{m π}{a}) s i n (\frac{n π}{b})}{ω_{m n}} [\begin{array}{l} (s i n (\frac{m π x}{a}) c o s (\frac{m π}{2}) + c o s (\frac{m π x}{a}) s i n (\frac{m π}{2})) \\ \times (s i n (\frac{n π}{b}) c o s (\frac{n π}{2}) + c o s (\frac{n π y}{a}) s i n (\frac{n π}{2})) \end{array}]$ (57)

The exact analytical solution for the natural frequency is given as

$ω_{m n} = \sqrt{\frac{D_{r}}{ρ h}} [{(\frac{m π}{a})}^{2} + {(\frac{n π}{b})}^{2}]$ $ω_{m n} = \sqrt{\frac{D_{r}}{ρ h}} [{(\frac{m π}{a})}^{2} + {(\frac{n π}{b})}^{2}]$ (58)

Parameter used for simulations

Case 1:

Plate dimension: 1m x 1m

Plate thickness: 0.025 m, 0.05 m and 0.075 m

Point load: 4500 N

Uniformly distributed load = 4.5 kN/m

The mechanical properties of the plate were E= 205GPa and ρ= 8050 kg/m3. The Poisson’s ratio is 0.3 as shown in Table 1a.

Material	Density kg/m³	Yield stress (Mpa)	Young's modulus (N/m2)	Poisson's ratio
Steel	8050	250	205x10⁹	0.3

Table 1a Physical properties of the materials used for the simulationsp

Case 2:

Plate dimension: 700 mm x 700 mm

Plate thickness: 10 mm

Point load: 800 N

The mechanical properties of the plate were E= 205GPa and ρ= 8050 kg/m3. The Poisson’s ratio is 0.3 as shown in Table 1b.

Material	Density kg/m³	Yield stress (Mpa)	Young's modulus (N/m²)	Poisson's ratio
Aluminium	2700	276	70x10⁹	0.3

Table 1b Physical properties of the materials used for the simulations

Result and discussion

The results of the numerical simulations using MATLAB and FlexPDE are presented in this section. Grid independence and sensitivity analyses are carried out. Figures 6–8 present the finite element discretization of the plate.

Figure 6 Discretisation of 1 m by 1 m Mild steel plate.

Figure 7 Triangular element mesh and refined mesh (MATLAB) for a rectangular plate of length = 500 mm, width = 200 mm and thickness = 12 mm.

Figure 8 Triangular element mesh and refined mesh (MATLAB) for a rectangular plate of length = 500 mm, width = 300 mm and thickness = 12 mm.

Grid independency test and code verification

A mesh sensitivity analysis was carried out to ensure grid independence. In order to choose the grid size, grid independency test is performed for different grid mesh of sizes. The number of grid size from which there is no appreciable change in the results and the accuracy is not affected even if the numbers of elements are increased by decreasing the size of the elements.

In order to verify and establish the reliability of the finite element methods used in this work, we made some comparisons of the results of an in-house MATLAB code and that of commercial computational solid mechanics code using FLEXPDE as well as exact analytical method. It is established that, the two finite element methods have been proven to be adequate and give accurate results for vibration problem presented in this project. The results show excellent agreement between the exact solution and the finite element solution.

Clamped-clamped-clamped-simply supported (CCCS) plate subjected to point load

Table 2 shows the effect of plate thickness on the deflection of the plate for the Clamped-Clamped-Clamped-Simply Supported (CCCS) boundary conditions subjected to point load of 4.5kN.

	Maximum deflection of the plate (m)
S. No. Thickness	Exact analytical solution	FEM MATLAB	FEM FLEXPDE
S. No. Thickness		FEM MATLAB	FEM FLEXPDE
1 0.025	4.19×10^-2	4.86×10^-2	5.45×10^-2
2 0.050	6.58×10^-3	7.45×10^-3	7.93×10^-3
3 0.075	2.23×10^-3	3.04×10^-3	2.79×10^-3

Table 2 Comparison of the deflection results for plate of length = 1 m, width = 1 m and Point load = 4500 N

Clamped-clamped-clamped-simply supported (CCCS) plate subjected to uniformly distributed load

While Table 3 shows the effect of plate thickness on the deflection of the plate for the same boundary conditions subjected to a uniformly distributed load of 4.5kN/m.

	Maximum deflection of the plate (m)
S. No. Thickness	Exact analytical solution	FEM MATLAB	FEM FLEXPDE
S. No. Thickness		FEM MATLAB	FEM FLEXPDE
1 0.025	5.65×10^-3	5.86×10^-3	5.79×10^-3
2 0.050	7.98×10^-4	8.45×10^-4	8.20×10^-4
3 0.075	2.53×10^-4	3.04×10^-4	2.79×10^-4

Table 3 Comparison of the deflection results for plate of length = 1 m, width = 1 m and uniformly distributed load = 4.5kN/m²

The results obtained for deflections from FLEXPDE give good agreement with the results obtained from FEM MATLAB with the FLEXPDE values slightly higher. However, the difference between the two values reduces from at thickness of the plate increases. Also, it was observed that the deflection of the plate decreases with increase in the thickness of the plate. This is because no shear deformation is considered in the analysis of thin plate. If the plate is thick the effect of shear deformation becomes significant and therefore the result obtained may vary. Therefore, to reduce deflection, the thickness of the plate should be increased.

Clamped-simply supported-clamped-simply supported (CSCS) plate subjected to point load

Table 4 shows the effect of plate thickness on the deflection of the plate for the Clamped- Clamped-Simply Supported -Clamped-Simply Supported (CSCS) boundary conditions subjected to point load of 4.5kN.

	Maximum deflection of the plate (m)
S. No. Thickness	Exact analytical solution	FEM MATLAB	FEM FLEXPDE
S. No. Thickness		FEM MATLAB	FEM FLEXPDE
1 0.025	5.23×10^-2	6.17×10^-2	5.45×10^-2
2 0.050	7.34×10^-3	8.23×10^-3	8.52×10^-3
3 0.075	2.51×10^-3	2.78×10^-3	2.83×10^-3

Table 4 Comparison of the deflection results for plate of length = 1 m, width = 1 m and Point load = 4500 N

Clamped-simply supported-clamped-simply supported (CSCS) plate subjected to uniformly distributed load

While Table 5 shows the effect of plate thickness on the deflection of the plate for the same boundary conditions subjected to a uniformly distributed load of 4.5kN/m.

	Maximum deflection of the plate (m)
S. No. Thickness	Exact analytical solution	FEM MATLAB	FEM FLEXPDE
S. No. Thickness		FEM MATLAB	FEM FLEXPDE
1 0.025	7.65×10^-3	7.55×10^-3	8.69×10^-3
2 0.050	1.22×10^-3	1.12×10^-3	1.15×10^-3
3 0.075	3.46×10^-4	3.44×10^-4	3.56×10^-4

Table 5 Comparison of the deflection results for plate of length = 1 m, width = 1 m and uniformly distributed load = 5kN/m²

Clamped -clamped- simply supported -simply supported (CCSS) plate subjected to point load

Table 6 shows the effect of plate thickness on the deflection of the plate for the Clamped- Clamped-Clamped-Simply Supported-Simply Supported (CCSS) boundary conditions subjected to point load of 4.5kN.

	Maximum deflection of the plate (m)
S. No. Thickness	Exact analytical solution	FEM MATLAB	FEM FLEXPDE
S. No. Thickness		FEM MATLAB	FEM FLEXPDE
1 0.025	5.67×10^-2	6.17×10^-2	5.45×10^-2
2 0.050	8.47×10^-3	8.51×10^-3	8.58×10^-3
3 0.075	2.71×10^-3	2.82×10^-3	2.86×10^-3

Table 6 Comparison of the deflection results for plate of length = 1 m, width = 1 m and Point load = 4500 N

Clamped -clamped-simply supported -simply supported (CCSS) plate subjected to uniformly distributed load

While Table 7 shows the effect of plate thickness on the deflection of the plate for the same boundary conditions subjected to a uniformly distributed load of 4.5kN/m.

	Maximum deflection of the plate (m)
S. No. Thickness	Exact analytical solution	FEM MATLAB	FEM FLEXPDE
S. No. Thickness		FEM MATLAB	FEM FLEXPDE
1 0.025	8.27×10^-3	9.55×10^-3	9.71×10^-3
2 0.050	1.15×10^-3	1.24×10^-3	1.26×10^-3
3 0.075	3.52×10^-4	3.61×10^-4	3.70×10^-4

Table 7 Comparison of the deflection results for plate of length = 1 m, width = 1 m and uniformly distributed load = 5kN/m²

Clamped -simply supported- simply supported -simply supported (CSSS) plate subjected to point load

Tables 8 shows the effect of plate thickness on the deflection of the plate for the Clamped -Simply Supported-Simply Supported -Simply Supported (CCSS) boundary conditions subjected to point load of 4.5kN.

	Maximum deflection of the plate (m)
S. No. Thickness	Exact analytical solution	FEM MATLAB	FEM FLEXPDE
S. No. Thickness		FEM MATLAB	FEM FLEXPDE
1 0.025	7.25×10^-2	6.89×10^-2	6.85×10^-2
2 0.050	1.02×10^-2	9.47×10^-3	9.49×10^-3
3 0.075	3.15×10^-3	3.13×10^-3	3.14×10^-3

Table 8 Comparison of the deflection results for plate of length = 1 m, width = 1 m and Point load = 4500 N

Clamped -simply supported- simply supported -simply supported (CSSS) plate subjected to uniformly distributed load

While Table 9 shows the effect of plate thickness on the deflection of the plate for the same boundary conditions subjected to a uniformly distributed load of 4.5kN/m.

	Maximum deflection of the plate (m)
S. No. Thickness	Exact analytical solution	FEM MATLAB	FEM FLEXPDE
S. No. Thickness		FEM MATLAB	FEM FLEXPDE
1 0.025	1.07×10^-2	1.15×10^-3	1.15×10^-3
2 0.050	1.45×10^-3	1.50×10^-3	1.50×10^-3
3 0.075	4.31×10^-4	4.71×10^-4	4.70×10^-4

Table 9 Comparison of the deflection results for plate of length = 1 m, width = 1 m and uniformly distributed load = 5kN/m²

Figures 9&10 shows the effects of plate thickness and boundary conditions on the deflection of the plate under point load. While Figures 11–16 presents the impacts of the plate thickness and boundary conditions on the deflection of the plate under uniformly distributed load. The figures show that the deflection of the plate decreases as the thickness of the plate increases. However, the higher the number of simply supported boundaries, the higher the deflection of the plate. This also connotes that the higher the number of clamped-supported boundaries, the lower the deflection of the plate. Figures 17–21 shows the three-dimensional view of the vibration of the plate using FLEXPDE.

Figure 9 Deflection versus dimension plot of point distributed load under CCCS boundary conditions.

Figure 10 Deflection versus dimension plot of point distributed load under CSCS boundary conditions.

Figure 11 Deflection versus dimension plot of uniformly distributed load under CCSS boundary conditions.

Figure 12 Deflection versus dimension plot of uniformly distributed load under CSSS boundary conditions.

Figure 13 Deflection versus dimension plot of uniformly distributed load under CCCS boundary conditions.

Figure 14 Deflection versus dimension plot of uniformly distributed load under CSCS boundary conditions.

Figure 15 Deflection versus dimension plot of uniformly distributed load under CCSS boundary conditions.

Figure 16 Deflection versus dimension plot of uniformly distributed load under CSSS boundary conditions.

Figure 17 Different modes of vibration (FLEXPDE FEM) for a rectangular plate of length = 1000 mm, width = 1000 mm and thickness = 75 mm.

Figure 18 Different modes of vibration (FLEXPDE FEM) for a rectangular plate of length = 200 mm, width = 100 mm and thickness = 12 mm.

Figure 19 Different modes of vibration (FLEXPDE FEM) for a rectangular plate of length = 500 mm, width = 100 mm and thickness = 12 mm.

Figure 20 Different modes of vibration (FLEXPDE FEM) for a rectangular plate of length = 500 mm, width = 100 mm and thickness = 12 mm.

Figure 21 Different modes of vibration (FLEXPDE FEM) for a rectangular plate of length = 700 mm, width = 700 mm and thickness = 10 mm.

The natural frequencies of plates are investigated versus modal numbers by varying the length and width of the plates with clamped-clamped-clamped - simply supported (CCCS), clamped-simply-clamped - simply supported (CSCS), clamped-clamped-simply - simply supported (CCSS), clamped-simply-simply- simply supported (CSSS) boundary conditions as presented in Table 10. The frequencies of the plates increase by increasing the modal number, and CCCS frequencies are greater than the frequencies of other boundary conditions. It is observed from the Table that the boundary conditions have a significant influence on the vibration of plates. Also, the obtained results clearly show that frequency parameters increase if more constraints are included. Simply supported at all edges have lower frequency than clamped edge due to clamping in all edges. This means that as the constraint on the edges increases the flexural rigidity of the plate increases and hence there is an increase in the frequency. Table 11 presents comparison of results of the first 6 natural frequencies and mode shapes of simply supported rectangular plates.

Modes	CCCS	CSCS	CCSS	CSSS
(1, 1)	54.34	38.45	41.84	10.56
(1, 2)	55.13	40.11	42.48	15.98
(2, 2)	56.06	42.97	44.25	27.53
(1, 3)	64.98	53.65	55.31	32.87
(2, 3)	68.75	57.85	59.02	39.21
(1, 4)	75.27	67.33	69.18	48.9
(3, 3)	82.81	74.72	76.89	57.42

Table 10 Effects of boundary conditions on the natural frequencies at different mode numbers

Mode number	6.25 mm	12.5 mm	25 mm	50 mm
1	136.9	272.2	535.3	1024
2	263.3	522.27	1019.8	1906.8
3	422.3	833.6	1634.9	3008.8
4	474.5	940.3	1823.1	3100.6
5	547.9	1083.3	2088.3	3316
6	758.1	1493.7	2847.7	3748.1

Table 11 Effects of plate thickness conditions on the natural frequencies at different mode numbers (SSSS)

The table shows that the results of the natural frequency in this work agree very well with the results of the natural frequencies of the exact analytical solution and those numerical solutions available in literature and with exact solutions. Also, it was established that accuracy in the finite element solution increases as the thickness of the rectangular plate increases.

Conclusion

In this work, finite element method is applied to structural analysis of rectangular plate subjected to different boundary conditions. The verification and the reliability of the finite element methods used in this work have been established by comparing the results of an in-house MATLAB code and that of commercial computational solid mechanics code using FLEXPDE as well as exact analytical method. The results of the exact analytical method, MATLAB finite element method and FLEX finite element method show excellent agreements. Also, it was observed that the boundary conditions and the dimensions of the plate have significant influence on the vibration of plates. The results showed that Clamped- Simply Supported - Simply Supported -Simply Supported (CSSS) mixed boundary condition gave the highest deflection for each thickness and loading type while the Clamped-Clamped- Clamped-Simply Supported (CCCS) mixed boundary condition provided the lowest deflection. The FEM algorithm developed in this project shows great agreement with commercial software results and can be used by any person who has the MATLAB code. Although, the processing time for MATLAB symbolic language is higher than if one uses numbers in the equations. However, the authors believe that using the first strategy is better for code learning as the code consists of a good way to learn the finite element method in its state-of-art. It could be inferred from the study that use of the developed finite element codes in MATLAB is an easier, cheaper and quicker approach for the analysis of free vibration of a thin isotropic rectangular plate.

Acknowledgments

We thank University of Lagos, Nigeria and Lagos State University, Nigeria for providing us conducive environment and essential materials to carry out this research.