Investigation of mechanical responses of flexo-magnetic variable thickness nanoplates resting on elastic foundations, taking into account geometrical imperfections

Abstract

This work utilizes the finite element approach together with an innovative shear strain theory to investigate the static bending behavior, free vibration features, and static buckling phenomena of flexo-magnetic nanoplates. The inquiry specifically examines the fluctuation in both the thickness of the plate and the elasticity of the foundation. The influence of initial geometrical imperfections, including several categories such as local and global faults, is also taken into account. The influences of several factors, including the law governing thickness fluctuation, types of imperfections, boundary conditions, and elastic foundation, on the mechanical response of the plate are considered. Outcomes of the work include new and original discoveries that have not been discussed in previous research, adding to both theoretical comprehension and practical implementation.

Graphical abstract

Keywords

bending / vibration / buckling / nanoplates / imperfect / foundation

1 Introduction

Nanostructures have become more common in several industries and the flexomagnetic effect has garnered considerable interest, being applied in several kinds of electrical components. The mechanical behavior of these structures is of great importance in practical applications, especially in aiding the logical design of devices. In recent decades, a multitude of novel materials and intelligent gadgets have emerged. These include piezoelectric materials, such as those exhibiting flexoelectric and flexomagnetic effects. These materials have the capacity to store energy and can use mechanical deformation as a source of energy, reversibly, such as for electrical and magnetic effects [1]. This technology is very compatible with compact electronic devices, such as sensors and components used in electronic circuits. Scholars have conducted research on the mechanical responses of nanostructures composed of materials exhibiting electromechanical or electrochemical interactions. Liu et al. [2] focused on nonlocal and Kirchhoff theories for piezoelectric nanoplate thermo-electro-mechanical free oscillation. In that work, the piezoelectric nanoplate was considered to be a rectangular plate supported by a biaxial force, with an external electric voltage, and homogeneous thermal change. Hamilton’s principle was used to construct the governing equations and boundary constraints, which were solved analytically to determine the piezoelectric nanoplate’s natural frequencies. Zhang et al. [3] investigated the dispersion properties of elastic waves that propagate through a mono-layer nanoplate, taking into account both the surface piezoelectricity and the nonlocal small-scale effect. The general governing equations were derived using nonlocal electroelasticity theory, which involved the introduction of an intrinsic length. Additionally, the boundary conditions of the piezoelectric nanoplate were considered by incorporating the surface piezoelectric model and the generalized Young-Laplace equations to account for surface effects. Ke et al. [4] studied the thermo-electro-mechanical oscillation of a rectangular piezoelectric nanoplate, considering different boundary constraints. The investigation was based on the nonlocal hypothesis and the Mindlin plate theory. Most recently, Li et al. [5] examined the size-dependent wave propagation characteristics under thermoelectric loading of porous functionally graded piezoelectric (FGP) nanoplates deposited on a viscoelastic substrate. Van Minh and Van Ke [6] conducted an extensive investigation into the static bending, free oscillation, and dynamic responses of piezoelectric nanoplates with non-uniform thicknesses. These nanoplates were placed on elastic foundations that varied in their properties, and the flexoelectric effect was considered. Additionally, the thickness of the nanoplates varied in both the length and breadth directions. Zhang et al. [7] presented properties of size-dependent dispersion in piezoelectric nanoplates involving the surface effect. Some of the most relevant works on this issue can be seen in the following articles [8–14].

In the course of panel structure production, the occurrence of protrusions leading to surface irregularities is a common phenomenon, often referred to as an imperfection. The presence of such flaws may significantly impact structural functionality, whereas investigations into the mechanical behavior of materials exhibiting faults in their basic form have yielded certain findings. Tho et al. [15] initially integrated the finite element method (FEM) with a newly developed third-order shear deformation beam theory to model the static bending and free oscillation behaviors of rotating piezoelectric nanobeams. The nanobeams were subjected to geometrical imperfections and flexoelectric effects while being supported by elastic foundations, based on Pasternak’s model. Li et al. [16] conducted a study wherein they disclosed and juxtaposed the impacts of the nonlocal scale, poor interface, interlaminar van der Waals force, and loading ratio. The imperfect interface was characterized using normal and shear springs, while the van der Waals force was quantified using the Hamaker formula. Sahmani and Aghdam [17] investigated the temperature-dependent nonlocal instability of hybrid functionally graded material exponential shear deformable nanoshells, considering their susceptibility to imperfections. Alam and Mishra [18] explored the postcritical imperfection sensitivity of FGP cylindrical nanoshells employing a boundary layer solution. The use of plates of varying thickness in some devices, either to optimize size or owing to spatial constraints, has prompted scientific investigation into the behavior of such plates. Duc et al. [19] analyzed the oscillation behavior and static buckling behavior of variable flexoelectric nanoplates. This investigation used the FEM in conjunction with a unique shear deformation theory based on hyperbolic sine functions. The thickness of the nanoplates was manipulated using both linear and nonlinear rules. Van Minh and Van Ke [6] considered the influence of plate thickness with different changing rules. Wang et al. [20] examined surface stress size dependence in nonlinear free vibration of FG quasi-3D nanoplates with arbitrary geometries and varying thicknesses using isogeometric analysis (IGA). Related works can be found in references [21–35]. Refs. [36–40] recently used the isogeometric analysis method to compute and enhance the performance of piezoelectric and flexoelectric materials.

This study primarily investigates the static bending, buckling, and oscillation reactions of nanoplates with varying thicknesses, considering the influence of the flexomagnetic effect. Additionally, the nanoplates are supported by an elastic foundation, with the foundation parameter also being subject to change. This research thus addresses a novel aspect that has not been previously explored, presenting several intriguing aspects that are not discussed in existing literature. The primary contribution of this study is in the provision of calculation formulae and comprehensive surveys, which serve as a foundation for the design and use of nanoplate structures with optimal efficiency. This work aims to enhance the practicality and utility of the research in this field.

The present work is organized into five different sections. Section 2 presents calculation methods for three distinct issues, namely static bending, specific vibration, and buckling of nanoplates, with changing thickness. Section 3 presents illustrative examples that serve to demonstrate the dependability of the calculation theory outlined in Section 2. Additionally, this section includes the presentation and discussion of numerical calculation results. The findings of the investigation are presented in Section 4.

2 Calculation formulations for nanoplates with variable thickness, taking into account initial geometrical imperfections

Consider a variable-thickness nanoplate with length and width parameters a and b, and the thickness varies in terms of both x and y coordinates as shown in Fig.1.

In this study, an enhanced shear deformation theory is used to derive the calculation formula for plates. Therefore, the displacement field along the three coordinate axes Ox, Oy, and Oz are denoted as v_x, v_y, and v_z accordingly, and may be represented as follows [41]:

(1)\{\begin{array}{l}{v}_x=-\mathit{\text{z}}{v}_{\mathit{\text{z}}\mathrm{b},x}-f(\mathit{\text{z}})v_{\mathit{\text{z}}\mathrm{s},x},\\ v_y=-\mathit{\text{z}}{v}_{\mathit{\text{z}}\mathrm{b},y}-f(\mathit{\text{z}})v_{\mathit{\text{z}}\mathrm{s},y},\\ v_{\mathit{\text{z}}}=v_{\mathit{\text{z}}\mathrm{b}}+v_{\mathit{\text{z}}\mathrm{s}},\end{array}

where the comma at the subscript denotes the derivative with respect to the variable immediately following it; v_z is divided into two parts: v_zb and v_zs; and f(\mathit{\text{z}})=\mathit{\text{z}}- t(\mathit{\text{z}}), t(\mathit{\text{z}})=h\cdot \text{sinh}{\displaystyle \frac{\mathit{\text{z}}}{h}}-\mathit{\text{z}}\cdot \text{cosh}{\displaystyle \frac{1}{2}}.

Since initial shape imperfections are taken into account, variable components are calculated as:

(2)\{\begin{array}{l}\text{ε}={\left\{{\epsilon }_x,{\epsilon }_y,{\epsilon }_{xy}\right\}}^{\text{T}}=\mathit{\text{z}}{\text{ε}}_{\mathrm{b}}+r(\mathit{\text{z}}){\text{ε}}_{\mathrm{s}}+{\text{ε}}_{\text{im}},\\ \text{γ}={\left\{{\gamma }_{x\mathit{\text{z}}},{\gamma }_{y\mathit{\text{z}}}\right\}}^{\text{T}}=t_{,\mathit{\text{z}}}\left(\mathit{\text{z}}\right){\text{γ}}_0,\end{array}

where

(3)\begin{array}{rl}& {\epsilon }_x=-\mathit{\text{z}}{v}_{\mathit{\text{z}}\mathrm{b},xx}-r(\mathit{\text{z}})v_{\mathit{\text{z}}\mathrm{s},xx}=\mathit{\text{z}}{\epsilon }_{\mathrm{b}x}+r(\mathit{\text{z}}){\epsilon }_{\mathrm{s}x},\\ & {\epsilon }_y=-\mathit{\text{z}}{v}_{\mathit{\text{z}}\mathrm{b},yy}-r(\mathit{\text{z}})v_{\mathit{\text{z}}\mathrm{s},yy}=\mathit{\text{z}}{\epsilon }_{\mathrm{b}y}+r(\mathit{\text{z}}){\epsilon }_{\mathrm{s}y},\\ & {\epsilon }_{\mathrm{i}\mathrm{m}x}=v_{\mathit{\text{z}},x}{v}_{\mathrm{i}\mathrm{m},x},\\ & {\epsilon }_{\mathrm{i}\mathrm{m}y}=v_{\mathit{\text{z}},y}{v}_{\mathrm{i}\mathrm{m},y},\\ & {\gamma }_{\mathrm{i}\mathrm{m}xy}=v_{\mathit{\text{z}},x}{v}_{\mathrm{i}\mathrm{m},y}+v_{\mathit{\text{z}},y}{v}_{\mathrm{i}\mathrm{m},x},\\ & {\gamma }_{xy}=-2\mathit{\text{z}}{v}_{\mathit{\text{z}}\mathrm{b},xy}-2r(\mathit{\text{z}})v_{\mathit{\text{z}}\mathrm{s},xy}=\mathit{\text{z}}{\gamma }_{\mathrm{b}xy}^{}+r(\mathit{\text{z}}){\gamma }_{\mathrm{s}xy},\\ & {\gamma }_{x\mathit{\text{z}}}=t{\left(\mathit{\text{z}}\right)}_{,\mathit{\text{z}}}{v}_{\mathit{\text{z}}\mathrm{s},x},\\ & {\gamma }_{y\mathit{\text{z}}}=t{\left(\mathit{\text{z}}\right)}_{,\mathit{\text{z}}}{v}_{\mathit{\text{z}}\mathrm{s},y},\end{array}

where v_im is the original shape imperfection, it transforms in the Oxy plane according to Eq. (4):

(4)\begin{array}{rl}{v}_{\mathrm{i}\mathrm{m}}=& {\eta }_{\mathrm{i}\mathrm{m}}h\mathrm{s}\mathrm{e}\mathrm{c}\mathrm{h}\left[{\delta }_1\left(\frac{x}{a}-{\psi }_1\right)\right]\mathrm{c}\mathrm{o}\mathrm{s}\left[{\mu }_1\left(\frac{x}{a}-{\psi }_1\right)\right]\\ & \cdot \mathrm{s}\mathrm{e}\mathrm{c}\mathrm{h}\left[{\delta }_2\left(\frac{y}{b}-{\psi }_2\right)\right]\mathrm{c}\mathrm{o}\mathrm{s}\left[{\mu }_2\left(\frac{y}{b}-{\psi }_2\right)\right],\end{array}

where {\eta }_{\mathrm{i}\mathrm{m}} is the magnitude of the imperfection; {\delta }_1,{\psi }_1, {\mu }_1, {\delta }_2,{\psi }_2, and {\mu }_2 are the parameters representing the global or local influence of the imperfection. Depending on the parameters of the v_im function, the strain expression is different, so the mechanical response of the nanoplate is also different. Fig.2 shows some imperfect forms of the original shape, which are also the forms that will be used for calculations in this work.

The calculation of strain gradients involves the computation of individual derivatives:

(5)\begin{array}{rl}\mathit{{\rm Y}}=& \left\{\begin{array}{c}{\mathrm{{\rm Y}}}_{113}=-v_{\mathit{\text{z}}\mathrm{b},xx}-f_{,\mathit{\text{z}}}{v}_{\mathit{\text{z}}\mathrm{s},xx}\\ {\mathrm{{\rm Y}}}_{223}=-v_{\mathit{\text{z}}\mathrm{b},yy}-f_{,\mathit{\text{z}}}{v}_{\mathit{\text{z}}\mathrm{s},yy}\end{array}\right\}\\ =& \left\{\begin{array}{c}-v_{\mathit{\text{z}}\mathrm{b},xx}\\ -v_{\mathit{\text{z}}\mathrm{b},yy}\end{array}\right\}+f_{,\mathit{\text{z}}}\left(\mathit{\text{z}}\right)\left\{\begin{array}{c}-v_{\mathit{\text{z}}\mathrm{s},xx}\\ -v_{\mathit{\text{z}}\mathrm{s},yy}\end{array}\right\}={\mathit{{\rm Y}}}_{\mathrm{b}}+f_{,\mathit{\text{z}}}{\mathit{{\rm Y}}}_{\mathrm{s}}.\end{array}

It should be noted that in the general case, there is a strain gradient component along the Oz axis. However, this component is smaller than the two strain gradient components deforming in the Ox and Oy directions, so it is ignored [42].

The normal stress and shear stress components are expressed as follows [43–45]:

\begin{array}{rl}\text{σ}=& \left\{\begin{array}{c}{\sigma }_x\\ {\sigma }_y\\ {\tau }_{xy}\end{array}\right\}=\underset{{\mathit{C}}_{\epsilon }}{\underbrace{\left[\begin{array}{ccc}{c}_{11}& c_{12}& 0\\ c_{12}& c_{11}& 0\\ 0& 0& c_{66}\end{array}\right]}}\left\{\begin{array}{c}{\epsilon }_x\\ {\epsilon }_y\\ {\gamma }_{xy}\end{array}\right\}-\underset{{\mathit{Q}}_{31}}{\underbrace{{\overline{q}}_{31}\left\{\begin{array}{c}1\\ 1\\ 0\end{array}\right\}}}{M}_{\mathit{\text{z}}}\\ =& {\mathit{C}}_{\epsilon }\text{ε}-{\mathit{Q}}_{31}{\mathrm{M}}_{\mathit{\text{z}}},\end{array}

(6)\text{τ}=\left\{\begin{array}{l}{\tau }_{xz}\\ {\tau }_{yz}\end{array}\right\}=\left[\begin{array}{cc}{c}_{66}& \text{0}\\ 0& c_{66}\end{array}\right] \text{γ} ={\mathit{C}}_{\gamma }\text{γ},

where c_ij are material constants, {\overline{q}}_{31} is the piezomagnetic parameter, M_z is magnetic field.

The higher-order stress components are represented as follows [45]:

(7)\text{Ψ}=\left\{\begin{array}{l}{\Psi }_{xx\mathit{\text{z}}}\\ {\Psi }_{yy\mathit{\text{z}}}\end{array}\right\}=\underset{{\text{Θ}}_{14}}{\underbrace{\left[\begin{array}{cc}{\mathit{\text{ƛ}}}_{14}& 0\\ 0& {\mathit{\text{ƛ}}}_{14}\end{array}\right]}}\left\{\begin{array}{l}{{\rm Y}}_{xxz}\\ {{\rm Y}}_{yyz}\end{array}\right\} -\underset{{\widetilde{\text{Θ}}}_{14}}{\underbrace{f_{14}\left\{\begin{array}{l}1\\ 1\end{array}\right\}}} M_z= {\text{Θ}}_{14}{\rm Y}-{\widetilde{\text{Θ}}}_{14}{M}_{\mathit{\text{z}}} ,

where f₁₄ is the flexomagnetic coefficient, and {\mathit{\text{ƛ}}}_{14} is calculated from the material parameters as {\mathit{\text{ƛ}}}_{14} = {\mathit{\text{ƛ}}}_3 + 2{\mathit{\text{ƛ}}}_4:

(8)\begin{array}{rl}& {\mathit{\text{ƛ}}}_3=\frac{8}{3}{\mathit{\text{ƛ}}}_2+\frac{2}{5}{\mathit{\text{ƛ}}}_5,\\ & {\mathit{\text{ƛ}}}_4=\frac{\mu }{3}\left({\mathrm{\Gamma }}_1^2+6{\mathrm{\Gamma }}_2^2\right),\\ & {\mathit{\text{ƛ}}}_2=\frac{\mu }{30}\left(27{\mathrm{\Gamma }}_0^2-4{\mathrm{\Gamma }}_1^2-15{\mathrm{\Gamma }}_2^2\right),\\ & {\mathit{\text{ƛ}}}_5=\frac{\mu }{3}\left({\mathrm{\Gamma }}_1^2-3{\mathrm{\Gamma }}_2^2\right),\end{array}

where \mu is Lamé parameter, {\mathrm{\Gamma }}_i (i = 0,1,2) represents the material length scales.

The magnetic flux (MF_z) [44,45] is as follow:

(9)M{F}_{\mathit{\text{z}}}=d_{33}{M}_{\mathit{\text{z}}} +{\overline{q}}_{31} \left({\epsilon }_x+{\epsilon }_y\right)+f_{14}\left({{\rm Y}}_{xx\mathit{\text{z}}}+{{\rm Y}}_{yy\mathit{\text{z}}}\right) .

To find the magnetic field component M_z, one must find the magnetic potential, the relationship between these two components is as follow:

(10)M_{\mathit{\text{z}}}= -\frac{\mathrm{\partial }\psi }{\mathrm{\partial }\mathit{\text{z}}}.

This work does not consider the external magnetic field, the condition of the magnetic flux component is {\displaystyle \frac{\mathrm{\partial }M{F}_{\mathit{\text{z}}}}{\mathrm{\partial }\mathit{\text{z}}}}=0. Substituting Eq. (10) into Eq. (9), Eq. (11) is obtained as follow:

(11)\begin{array}{rl}\text{ψ}=& -\frac{{\overline{q}}_{31}}{d_{33}}\{\left(v_{\mathit{\text{z}}\text{b},xx}+v_{\mathit{\text{z}}\mathrm{b},yy}\right)\frac{{\mathit{\text{z}}}^2}{2}\\ & +\left(\frac{{\mathit{\text{z}}}^2}{2}-h^2\mathrm{c}\mathrm{o}\mathrm{s}\mathrm{h}\left(\frac{\mathit{\text{z}}}{h}\right)+\frac{{\mathit{\text{z}}}^2}{2}\mathrm{c}\mathrm{o}\mathrm{s}\mathrm{h}\left(\frac{1}{2}\right)\right)\left(v_{\mathit{\text{z}}\mathrm{s},xx}+v_{\mathit{\text{z}}\mathrm{s},yy}\right)\}\\ & -\frac{f_{14}}{d_{33}}{f}_{\mathit{\text{z}}}\left(v_{\mathit{\text{z}}\mathrm{s},xx}+v_{\mathit{\text{z}}\mathrm{s},yy}\right)+d_1\mathit{\text{z}}+d_0,\end{array}

where d₀ and d₁ are constants. To determine the expression of M_z like Eq. (10), it is necessary to determine the constant d₁, and the condition is as follow:

(12)\{\begin{array}{l}{M{F}_{\mathit{\text{z}}}|}_{\mathit{\text{z}}-h/2}=0\\ {M{F}_{\mathit{\text{z}}}|}_{\mathit{\text{z}}-h/2}=0\end{array}.

Combining Eqs. (9)–(12), one finds the constant d₁, from which it is easy to find the Eq. (13) for the magnetic field:

(13)\begin{array}{rl}{M}_{\mathit{\text{z}}}=& \frac{{\overline{q}}_{31}}{d_{33}}\left\{\left(v_{\mathit{\text{z}}\mathrm{b},xx}+v_{\mathit{\text{z}}\mathrm{b},yy}\right)\mathit{\text{z}}+f\left(\mathit{\text{z}}\right)\left(v_{\mathit{\text{z}}\mathrm{s},xx}+v_{\mathit{\text{z}}\mathrm{s},yy}\right)\right\}\\ & -\frac{{\overline{q}}_{31}}{d_{33}}\left(v_{\mathit{\text{z}},x}{v}_{\mathrm{i}\mathrm{m},x}+v_{\mathit{\text{z}},y}{v}_{\mathrm{i}\mathrm{m},y}\right)\\ & +\frac{f_{14}}{d_{33}}{f}_{,\mathit{\text{z}}}\left(\mathit{\text{z}}\right)\left(v_{\mathit{\text{z}}\mathrm{s},xx}+v_{\mathit{\text{z}}\mathrm{s},yy}\right)+\frac{f_{14}}{d_{33}}\left(v_{\mathit{\text{z}}\mathrm{b},xx}+v_{\mathit{\text{z}}\mathrm{b},yy}\right).\end{array}

To find the equilibrium equation for the plate, this paper uses Hamilton principle, according to which one the following Eq. (14) applies:

(14)\underset{t_0}{\overset{t_1}{\int }}\left(\delta E_{\mathrm{p}\mathrm{l}\mathrm{a}\mathrm{t}\mathrm{e}}+\delta E_{\mathrm{f}\mathrm{o}\mathrm{u}\mathrm{n}\mathrm{d}\mathrm{a}\mathrm{t}\mathrm{i}\mathrm{o}\mathrm{n}}-\delta E_{\text{kinetic}}-\delta E_{\text{load}}\right)\mathrm{d}t=0 ,

where E_plate, E_foundation, E_kinetic, and E_load are the deformation energy of the structure, the energy of the elastic foundation, the kinetic energy and the work of the external force, respectively. They are determined as follows:

\begin{array}{rl}\delta E_{\mathrm{p}\mathrm{l}\mathrm{a}\mathrm{t}\mathrm{e}}=& \underset{V}{\int }\left(\delta {\text{ε}}^{\mathrm{T}}\text{σ}+\delta {\text{γ}}^{\mathrm{T}}\text{τ}+\delta {\mathrm{{\rm Y}}}^{\mathrm{T}}\text{ψ}-\delta M{F}_{\mathit{\text{z}}}{M}_{\mathit{\text{z}}}\right)\mathrm{d}V\\ =& \underset{S}{\int }\underset{-\frac{h_{xy}}{2}}{\overset{\frac{h_{xy}}{2}}{\int }}\left(\delta {\text{ε}}^{\mathrm{T}}{\mathit{C}}_{\epsilon }\text{ε}+\delta {\text{γ}}^{\mathrm{T}}\text{τ}-\delta {\text{ε}}^{\mathrm{T}}{\mathit{Q}}_{31}{M}_{\mathit{\text{z}}}\mathrm{d}\mathit{\text{z}}\right)\mathrm{d}x\mathrm{d}y\\ & +\underset{S}{\int }\left(\underset{-\frac{h_{xy}}{2}}{\overset{\frac{h_{xy}}{2}}{\int }}\delta {\left({\mathit{{\rm Y}}}_{\mathrm{b}}+f_{,\mathit{\text{z}}}\left(\mathit{\text{z}}\right){\mathit{{\rm Y}}}_{\mathrm{s}}\right)}^{\mathrm{T}}{\text{Θ}}_{14}\left({\mathit{{\rm Y}}}_{\mathrm{b}}+f_{,\mathit{\text{z}}}\left(\mathit{\text{z}}\right){\mathit{{\rm Y}}}_{\mathrm{s}}\right)\mathrm{d}\mathit{\text{z}}\right)\mathrm{d}x\mathrm{d}y\end{array}

(15)\begin{array}{rl}& -\underset{S}{\int }\left(\underset{-\frac{h_{xy}}{2}}{\overset{\frac{h_{xy}}{2}}{\int }}\delta {\left({\mathit{{\rm Y}}}_{\mathrm{b}}+f_{,\mathit{\text{z}}}{\mathit{{\rm Y}}}_{\text{s}}\right)}^{\mathrm{T}}{\widetilde{\text{Θ}}}_{14}{M}_{\mathit{\text{z}}}\mathrm{d}\mathit{\text{z}}\right)\mathrm{d}S\\ & -\underset{S}{\int }\left(\underset{-\frac{h_{xy}}{2}}{\overset{\frac{h_{xy}}{2}}{\int }}{d}_{33}\delta M_{\mathit{\text{z}}}\cdot M_{\mathit{\text{z}}}\mathrm{d}\mathit{\text{z}}\right)\mathrm{d}x\mathrm{d}y\\ & -\underset{S}{\int }\left(\underset{-\frac{h_{xy}}{2}}{\overset{\frac{h_{xy}}{2}}{\int }}{\overline{q}}_{31}\delta \left({\epsilon }_x+{\epsilon }_y\right)M_{\mathit{\text{z}}}\mathrm{d}\mathit{\text{z}}\right)\mathrm{d}x\mathrm{d}y\\ & -\underset{S}{\int }\left(\underset{-\frac{h_{xy}}{2}}{\overset{\frac{h_{xy}}{2}}{\int }}{f}_{14}\delta \left({{\rm Y}}_{xx\mathit{\text{z}}}+{{\rm Y}}_{yy\mathit{\text{z}}}\right)M_{\mathit{\text{z}}}\mathrm{d}\mathit{\text{z}}\right)\mathrm{d}x\mathrm{d}y.\end{array}

The \delta E_{\mathrm{f}\mathrm{o}\mathrm{u}\mathrm{n}\mathrm{d}\mathrm{a}\mathrm{t}\mathrm{i}\mathrm{o}\mathrm{n}} is related to the elastic foundation as follows:

(16)\delta E_{\mathrm{f}\mathrm{o}\mathrm{u}\mathrm{n}\mathrm{d}\mathrm{a}\mathrm{t}\mathrm{i}\mathrm{o}\mathrm{n}}=\underset{S}{\int }\left(\begin{array}{c}{k}_{\mathrm{w}0}\left(1+{\beta }_x{\displaystyle \frac{x^2}{a^2}}\right)\left(1+{\beta }_y{\displaystyle \frac{y^2}{b^2}}\right)\delta v_{\mathit{\text{z}}}{v}_{\mathit{\text{z}}}\\ +k_{\mathrm{s}}\left(v_{\mathit{\text{z}},x}\delta \left(v_{\mathit{\text{z}},x}\right)+v_{\mathit{\text{z}},y}\delta \left(v_{\mathit{\text{z}},y}\right)\right)\end{array}\right)\mathrm{d}x\mathrm{d}y,

where k_w0 and k_s are parameters of the elastic foundation; {\beta }_x and {\beta }_y are parameters that represent the variation of the elastic foundation along the x and y axes.

Component \delta E_{\mathrm{l}\mathrm{o}\mathrm{a}\mathrm{d}} is calculated from the external force acting on the plate as follow:

(17)\delta E_{\mathrm{l}\mathrm{o}\mathrm{a}\mathrm{d}}=\underset{V}{\int }\delta {\text{ε}}_{\mathrm{i}\mathrm{m}\mathrm{p}\mathrm{r}\mathrm{e}\mathrm{s}\mathrm{s}}^{\mathrm{T}}{\sigma }_0{\text{ε}}_{\mathrm{i}\mathrm{m}\mathrm{p}\mathrm{r}\mathrm{e}\mathrm{s}\mathrm{s}}\mathrm{d}V+ \underset{S}{\int }\delta v_{\mathit{\text{z}}}{F}_{\mathit{\text{z}}}\mathrm{d}x\mathrm{d}y,

(18){\text{ε}}_{\mathrm{i}\mathrm{m}\mathrm{p}\mathrm{r}\mathrm{e}\mathrm{s}\mathrm{s}}=\left\{\begin{array}{c}{v}_{\mathit{\text{z}},x}\\ v_{\mathit{\text{z}},y}\end{array}\right\},{\sigma }_0= \left[\begin{array}{cc}{\sigma }_x^0& 0\\ 0& {\sigma }_y^0\end{array}\right]=F_0\left[\begin{array}{cc}{\displaystyle \frac{1}{h_{xy}}}& 0\\ 0& {\displaystyle \frac{1}{h_{xy}}}\end{array}\right] ,

where the compression force in the plate plane is F₀ and the force acting perpendicularly to the plate plane is F_z.

Component \delta E_{\mathrm{k}\mathrm{i}\mathrm{n}\mathrm{e}\mathrm{t}\mathrm{i}\mathrm{c}} is related to the kinetic energy as follow:

(19)\delta E_{\mathrm{k}\mathrm{i}\mathrm{n}\mathrm{e}\mathrm{t}\mathrm{i}\mathrm{c}} =\underset{S}{\int } \underset{-\frac{h_{xy}}{2}}{\overset{\frac{h_{xy}}{2}}{\int }}\left(\delta {\left\{v_{x,t},v_{y,t},v_{\mathit{\text{z}},t}\right\}}^{\mathrm{T}}\rho \left\{v_{x,t},v_{y,t},v_{\mathit{\text{z}},t}\right\}\mathrm{d}\mathit{\text{z}}\right)\mathrm{d}x\mathrm{d}y,

where t is time.

The present study used a four-node quadrilateral element. Each node is composed of six components, as shown below:

(20){\mathit{v}}_e=\sum _{i=1}^m{\left\{v_{\mathit{\text{z}}\mathrm{b}i},v_{\mathit{\text{z}}\mathrm{s}i},{\left(v_{\mathit{\text{z}}\mathrm{b},x}\right)}_i,{\left(v_{\mathit{\text{z}}\mathrm{s},x}\right)}_i,{\left(v_{\mathit{\text{z}}\mathrm{b},y}\right)}_i,{\left(v_{\mathit{\text{z}}\mathrm{s},y}\right)}_i\right\}}^{\mathrm{T}},

where m = 4, and the first two components in Eq. (20) are as follows:

(21)\{\begin{array}{c}{v}_{\mathit{\text{z}}\mathrm{b}}=\sum _{i=1}^m\left[V_i,V_{i+1},V_{i+2}\right]{\left[v_{\mathit{\text{z}}\mathrm{b}i},{\left(v_{\mathit{\text{z}}\mathrm{b},x}\right)}_i,{\left(v_{\mathit{\text{z}}\mathrm{b},y}\right)}_i\right]}^{\mathrm{T}}={\mathit{V}}_{\mathrm{b}}{\mathit{v}}_e,\\ v_{\mathit{\text{z}}\mathrm{s}}=\sum _{i=1}^m\left[V_i,V_{i+1},V_{i+2}\right]{\left[v_{\mathit{\text{z}}\mathrm{s}i},{\left(v_{\mathit{\text{z}}\mathrm{s},x}\right)}_i,{\left(v_{\mathit{\text{z}}\mathrm{s},y}\right)}_i\right]}^{\mathrm{T}}={\mathit{V}}_{\mathrm{s}}{\mathit{v}}_e,\end{array}

where Vi are Hermitian functions.

The last four components in (20) are interpolated as follows:

(22)\{\begin{array}{rl}& u_{\mathit{\text{z}}\mathrm{b},x}=\sum _{i=1}^m\left\{V_{i,x}{v}_{\mathit{\text{z}}\mathrm{b}i}+V_{i+1,x}{\left(u_{\mathit{\text{z}}\mathrm{b},i}\right)}_i+V_{i+2,x}{\left(u_{\mathit{\text{z}}\mathrm{b},y}\right)}_i\right\}={\mathit{V}}_{\mathrm{b},x}{\mathit{v}}_e,\\ & u_{\mathit{\text{z}}\mathrm{b},xx}={\mathit{V}}_{b,xx}{\mathit{v}}_e,\\ & u_{\mathit{\text{z}}b,y}=\sum _{i=1}^m\left\{V_{i,y}{v}_{\mathit{\text{z}}\mathrm{b}i}+V_{i+1,y}{\left(u_{\mathit{\text{z}}\mathrm{b},x}\right)}_i+V_{i+2,y}{\left(u_{\mathit{\text{z}}\mathrm{b},y}\right)}_i\right\}={\mathit{V}}_{\mathrm{b},y}{\mathit{v}}_e,\\ & u_{\mathit{\text{z}}\mathrm{b},yy}={\mathit{V}}_{\mathrm{b},yy}{\mathit{v}}_e,\\ & u_{\mathit{\text{z}}\mathrm{s},x}=\sum _{i=1}^m\left\{V_{i,x}{v}_{\mathit{\text{z}}\mathrm{s}i}+V_{i+1,x}{\left(u_{\mathit{\text{z}}\mathrm{s},x}\right)}_i+V_{i+2,x}{\left(u_{\mathit{\text{z}}\mathrm{s},y}\right)}_i\right\}={\mathit{V}}_{\mathrm{s},x}{\mathit{v}}_{e}v,\\ & u_{\mathit{\text{z}}\mathrm{s},xx}={\mathit{V}}_{\mathrm{s},xx}{\mathit{v}}_e,\\ & u_{\mathit{\text{z}}s,y}=\sum _{i=1}^m\left\{V_{i,y}{v}_{\mathit{\text{z}}\mathrm{s}i}+V_{i+1,y}{\left(u_{\mathit{\text{z}}\mathrm{s}},x\right)}_i+V_{i+2,y}{\left(u_{\mathit{\text{z}}\mathrm{s},y}\right)}_i\right\}={\mathit{V}}_{\mathrm{s},y}{\mathit{v}}_e,\\ & u_{\mathit{\text{z}}\mathrm{s},yy}={\mathit{V}}_{\mathrm{s},yy}{\mathit{v}}_e.\end{array}

Equations (3) and (5) are represented as element node displacement vectors as follows:

{\text{ε}}_{\text{b}}=-\left[\begin{array}{c}{\mathit{V}}_{\mathrm{b},xx}\\ {\mathit{V}}_{\mathrm{b},yy}\\ 2{\mathit{V}}_{\mathrm{b},xy}\end{array}\right]{\mathit{v}}_e={\tilde{\mathit{V}}}_{\mathrm{b}} {\mathit{v}}_e,

{\text{ε}}_{\mathrm{s}}=- \left[\begin{array}{c}{\mathit{V}}_{\mathrm{s},xx}\\ {\mathit{V}}_{\mathrm{s},yy}\\ 2{\mathit{V}}_{\mathrm{s},xy}\end{array}\right]{\mathit{v}}_e={\tilde{\mathit{V}}}_{\mathrm{s}}{\mathit{v}}_e,

(23){\text{ε}}_{\mathrm{i}\mathrm{m}} =\left[\begin{array}{c}\left({\mathit{V}}_{\mathrm{b},x}+{\mathit{V}}_{\mathrm{s},x}\right)v_{\mathrm{i}\mathrm{m}\mathrm{p},x}\\ \left({\mathit{V}}_{\mathrm{b},y}+{\mathit{V}}_{\mathrm{s},y}\right)v_{\mathrm{i}\mathrm{m}\mathrm{p},y}\\ \left({\mathit{V}}_{\mathrm{b},x}+{\mathit{V}}_{\mathrm{s},x}\right)v_{\mathrm{i}\mathrm{m}\mathrm{p},y}+\left({\mathit{V}}_{\mathrm{b},y}+{\mathit{V}}_{\mathrm{s},y}\right)v_{\mathrm{i}\mathrm{m}\mathrm{p},x}\end{array}\right]{\mathit{v}}_e={\stackrel{\mathbf{~}}{\mathit{V}}}_{\mathrm{i}\mathrm{m}}{\mathit{v}}_e.

{\text{γ}}_0 =\left[\begin{array}{c}{\mathit{V}}_{\mathrm{s},x}\\ {\mathit{V}}_{\mathrm{s},y}\end{array}\right]{\mathit{v}}_e={\stackrel{\mathbf{~}}{\mathit{V}}}_{\gamma } {\mathit{v}}_e,

(24)\begin{array}{c}{\mathit{{\rm Y}}}_{\mathrm{b}}=-\left[\begin{array}{c}{\mathit{V}}_{\mathrm{b},xx}\\ {\mathit{V}}_{\mathrm{b},yy}\end{array}\right]{\mathit{v}}_e={\stackrel{\mathbf{~}}{\mathit{V}}}_{\mathit{{\rm Y}}\mathrm{b}}{\mathit{v}}_e\text{,}\\ {\mathit{{\rm Y}}}_{\mathrm{s}}=-\left[\begin{array}{c}{\mathit{V}}_{\mathrm{s},xx}\\ {\mathit{V}}_{\mathrm{s},yy}\end{array}\right]{\mathit{v}}_e={\stackrel{\mathbf{~}}{\mathit{V}}}_{\mathit{{\rm Y}}\mathrm{s}}{\mathit{v}}_e.\end{array}

Substituting Eq. (23) and Eq. (24) into Eq. (15), Eq. (25) is obtained as follow:

(25)\begin{array}{rl}\delta E_{\mathrm{p}\mathrm{l}\mathrm{a}\mathrm{t}\mathrm{e}}=& \delta {\mathit{v}}_e^{\mathrm{T}}\underset{S_e}{\int }\underset{-\frac{h_{xy}}{2}}{\overset{\frac{h_{xy}}{2}}{\int }}\{\left(\mathit{\text{z}}{\stackrel{\mathbf{~}}{\mathit{V}}}_{\mathrm{b}}^{\mathrm{T}}+f{\stackrel{\mathbf{~}}{\mathit{V}}}_{\mathrm{s}}^{\mathrm{T}}+{\stackrel{\mathbf{~}}{\mathit{V}}}_{\mathrm{i}\mathrm{m}\mathrm{p}}^{\mathrm{T}}\right)\\ & \cdot {\mathit{C}}_{\epsilon }\left(\mathit{\text{z}}{\stackrel{\mathbf{~}}{\mathit{V}}}_b+f{\stackrel{\mathbf{~}}{\mathit{V}}}_{\mathrm{s}}+{\stackrel{\mathbf{~}}{\mathit{V}}}_{\mathrm{i}\mathrm{m}\mathrm{p}}\right){+\tilde{\mathit{V}}}_{\gamma }^{\mathrm{T}}t{\mathit{C}}_{\gamma }{\stackrel{\mathbf{~}}{\mathit{V}}}_{\gamma }\}\mathrm{d}\mathit{\text{z}}\mathrm{d}x\mathrm{d}y{\mathit{v}}_e\\ & -\delta {\mathit{v}}_e^{\mathrm{T}}{\mathit{Y}}_{\mathrm{b}}{\mathit{v}}_e-\delta {\mathit{v}}_e^{\mathrm{T}}{\mathit{Y}}_{\mathrm{s}}{\mathit{v}}_e-\delta {\mathit{v}}_e^{\mathrm{T}}{\mathit{Y}}_{1\mathrm{i}\mathrm{m}}{\mathit{v}}_e-\delta {\mathit{v}}_e^{\mathrm{T}}{\mathit{Y}}_{2\mathrm{i}\mathrm{m}}{\mathit{v}}_e\\ & +\delta {\mathit{v}}_e^{\mathrm{T}}{\mathit{Y}}_{1\mathrm{{\rm Y}}}{\mathit{v}}_e-\delta {\mathit{v}}_e^{\mathrm{T}}{\mathit{Y}}_{2\mathrm{{\rm Y}}}{\mathit{v}}_e-\delta {\mathit{v}}_e^{\mathrm{T}}{\mathit{Y}}_{3\mathrm{{\rm Y}}}{\mathit{v}}_e+\delta {\mathit{v}}_e^{\mathrm{T}}{\mathit{Y}}_{1\mathrm{{\rm Y}}\mathrm{m}}{\mathit{v}}_e\\ & -\delta {\mathit{v}}_e^{\mathrm{T}}{\mathit{Y}}_{4\mathrm{{\rm Y}}}{\mathit{v}}_e-\delta {\mathit{v}}_e^{\mathrm{T}}{\mathit{Y}}_{4\mathrm{{\rm Y}}}{\mathit{v}}_e-\delta {\mathit{v}}_e^{\mathrm{T}}{\mathit{Y}}_{5\mathrm{{\rm Y}}}{\mathit{v}}_e-\delta {\mathit{v}}_e^{\mathrm{T}}{\mathit{Y}}_{6\mathrm{{\rm Y}}}{\mathit{v}}_e\\ & +\delta {\mathit{v}}_e^{\mathrm{T}}{\mathit{Y}}_{2\mathrm{{\rm Y}}\mathrm{m}}{\mathit{v}}_e-\delta {\mathit{v}}_e^{\mathrm{T}}{\mathit{Y}}_{1\mathrm{b}\mathrm{s}}{\mathit{v}}_e-\delta {\mathit{v}}_e^{\mathrm{T}}{\mathit{Y}}_{2\mathrm{b}\mathrm{s}}{\mathit{v}}_e-\delta {\mathit{v}}_e^{\mathrm{T}}{\mathit{Y}}_{1\mathrm{b}\mathrm{s}\mathrm{i}\mathrm{m}}{\mathit{v}}_e\\ & +\delta {\mathit{v}}_e^{\mathrm{T}}{\mathit{Y}}_{2\mathrm{b}\mathrm{s}\mathrm{i}\mathrm{m}}{\mathit{v}}_e+\delta {\mathit{v}}_e^{\mathrm{T}}{\mathit{Y}}_{3\mathrm{b}\mathrm{s}\mathrm{i}\mathrm{m}}{\mathit{v}}_e+\delta {\mathit{v}}_e^{\mathrm{T}}{\mathit{Y}}_{3\mathrm{b}\mathrm{s}}{\mathit{v}}_e-\delta {\mathit{v}}_e^{\mathrm{T}}{\mathit{Y}}_{4\mathrm{b}\mathrm{s}\mathrm{i}\mathrm{m}}{\mathit{v}}_e\\ & -\delta {\mathit{v}}_e^{\mathrm{T}}{\mathit{Y}}_{5\mathrm{b}\mathrm{s}\mathrm{i}\mathrm{m}}{\mathit{v}}_e+\delta {\mathit{v}}_e^{\mathrm{T}}{\mathit{Y}}_{6\mathrm{b}\mathrm{s}\mathrm{i}\mathrm{m}}{\mathit{v}}_e-\delta {\mathit{v}}_e^{\mathrm{T}}{\mathit{Y}}_{7\mathrm{b}\mathrm{s}\mathrm{i}\mathrm{m}}{\mathit{v}}_e+\delta {\mathit{v}}_e^{\mathrm{T}}{\mathit{Y}}_{4\mathrm{b}\mathrm{s}}{\mathit{v}}_e\\ & -\delta {\mathit{v}}_e^{\mathrm{T}}{\mathit{Y}}_{8\mathrm{b}\mathrm{s}\mathrm{i}\mathrm{m}}{\mathit{v}}_e+\delta {\mathit{v}}_e^{\mathrm{T}}{\mathit{Y}}_{5\mathrm{b}\mathrm{s}}{\mathit{v}}_e=\delta {\mathit{v}}_e^{\mathrm{T}}\left({\mathit{K}}_e^{\mathrm{p}\mathrm{l}\mathrm{a}\mathrm{t}\mathrm{e}}\right){\mathit{v}}_e,\end{array}

where the coefficient matrices are shown as in Appendix in Electronic Supplementary Materials.

Equation (16) is rewritten as follow:

(26)\begin{array}{rl}& \delta E_{\mathrm{f}\mathrm{o}\mathrm{u}\mathrm{n}\mathrm{d}\mathrm{a}\mathrm{t}\mathrm{i}\mathrm{o}\mathrm{n}}\\ =& \delta {\mathit{v}}_e^{\mathrm{T}}\underset{S_e}{\int }(k_{\mathrm{w}0}\left(1+{\beta }_x\frac{x^2}{a^2}\right)\left(1+{\beta }_y\frac{y^2}{b^2}\right){\left({\mathit{V}}_{\mathrm{b}}+{\mathit{V}}_{\mathrm{s}}\right)}^{\mathrm{T}}\left({\mathit{V}}_{\mathrm{b}}+{\mathit{V}}_{\mathrm{s}}\right)\\ & +k_{\mathrm{s}}({\left({\mathit{V}}_{\mathrm{b},x}+{\mathit{V}}_{\mathrm{s},x}\right)}^{\mathrm{T}}\left({\mathit{V}}_{\mathrm{b},x}+{\mathit{V}}_{\mathrm{s},x}\right)\\ & +{\left({\mathit{V}}_{\mathrm{b},y}+{\mathit{V}}_{\mathrm{s},y}\right)}^{\mathrm{T}}\left({\text{V}}_{\mathrm{b},y}+{\mathit{V}}_{\mathrm{s},y}\right)))\mathrm{d}x\mathrm{d}y{\mathit{v}}_e\\ =& \delta {\mathit{v}}_e^{\mathrm{T}}{\mathit{K}}_e^{\mathrm{f}\mathrm{o}\mathrm{u}\mathrm{n}\mathrm{d}\mathrm{a}\mathrm{t}\mathrm{i}\mathrm{o}\mathrm{n}}{\mathit{v}}_e.\end{array}

Possible work Eq. (17) of external force is rewritten as follow:

(27)\begin{array}{rl}\delta E_{\mathrm{l}\mathrm{o}\mathrm{a}\mathrm{d}}=& \delta {\mathit{v}}_e^{\mathrm{T}}\underset{S_e}{\int }\left(\underset{-\frac{h_{xy}}{2}}{\overset{\frac{h_{xy}}{2}}{\int }}\left\{\left({\mathit{V}}_{\mathrm{b}}^{\mathrm{T}}+{\mathit{V}}_{\mathrm{s}}^{\mathrm{T}}\right){\sigma }_0\left({\mathit{V}}_{\mathrm{b}}+{\mathit{V}}_{\mathrm{s}}\right)\right\}\mathrm{d}\mathit{\text{z}}\right)\mathrm{d}x\mathrm{d}y{\mathit{v}}_e\\ & +\delta {\mathit{v}}_e^{\mathrm{T}}\underset{S_e}{\int }\left({\mathit{V}}_{\mathrm{b}}^{\mathrm{T}}+{\mathit{V}}_{\mathrm{s}}^{\mathrm{T}}\right)F_{\mathit{\text{z}}}\mathrm{d}x\mathrm{d}y=\delta {\mathit{v}}_e^{\mathrm{T}}{\mathit{K}}_{\mathrm{g}e}{\mathit{v}}_e^{}+\delta {\mathit{v}}_e^{\mathrm{T}}{\tilde{\mathit{F}}}_e.\end{array}

Equation (19) is written as follows:

(28)\delta E_{\mathrm{k}\mathrm{i}\mathrm{n}\mathrm{e}\mathrm{t}\mathrm{i}\mathrm{c}} =\delta {\mathit{v}}_{e,t}^{\mathrm{T}} \underset{S_e}{\int }\left(\underset{-\frac{h_{xy}}{2}}{\overset{\frac{h_{xy}}{2}}{\int }}\left({\mathit{V}}^{\mathrm{T}}{\mathit{G}}^{\mathrm{T}}\rho \mathit{G}\mathit{V}\right)\mathrm{d}\mathit{\text{z}}\right)\mathrm{d}x\mathrm{d}y{\mathit{v}}_{e,t}=\delta {\mathit{v}}_{e,t}^{\mathrm{T}} {\mathit{M}}_e{\mathit{v}}_{e,t},

where

(29)\mathit{G}=\left[\begin{array}{cccccc}0& 0& -\mathit{\text{z}}& -f(\mathit{\text{z}})& 0& 0\\ 0& 0& 0& 0& -\mathit{\text{z}}& -f(\mathit{\text{z}})\\ 1& 1& 0& 0& 0& 0\end{array}\right].

Combining Eq. (27) and Eq. (14), one obtains the balanced equation of the structure for the problem of static buckling as follows:

(30)\sum _e\left({\mathit{K}}_e^{\mathrm{p}\mathrm{l}\mathrm{a}\mathrm{t}\mathrm{e}}+{\mathit{K}}_e^{\mathrm{f}\mathrm{o}\mathrm{u}\mathrm{n}\mathrm{d}\mathrm{a}\mathrm{t}\mathrm{i}\mathrm{o}\mathrm{n}}-F_0{\mathit{K}}_{\mathrm{g}e}\right){\mathit{v}}_e=0 .

Consequently, the element stiffness matrix K_e and the geometrical stiffness matrix K_ge are influenced by the mechanical characteristics of the material and the plate’s geometry. Specifically, they are dependent on the law governing the thickness variation and the presence of form imperfections. These parameters have an impact on the plate’s critical buckling load as well.

The static balance equation of the plate is determined from Eqs. (25), (27), and (14) as follows:

(31)\sum _e\left({\mathit{K}}_e^{\mathrm{p}\mathrm{l}\mathrm{a}\mathrm{t}\mathrm{e}}+{\mathit{K}}_e^{\mathrm{f}\mathrm{o}\mathrm{u}\mathrm{n}\mathrm{d}\mathrm{a}\mathrm{t}\mathrm{i}\mathrm{o}\mathrm{n}}\right){\mathit{v}}_e=\sum _e{\tilde{\mathit{F}}}_e.

The dependence of the bending displacement of the plate on the change law of thickness and imperfections is attributed to the parameter known as the stiffness matrix K_e.

As for the problem of specific vibrations of the plate, the equation has the form as follow:

(32)\sum _e\left({\mathit{K}}_e^{\mathrm{p}\mathrm{l}\mathrm{a}\mathrm{t}\mathrm{e}}+{\mathit{K}}_e^{\mathrm{f}\mathrm{o}\mathrm{u}\mathrm{n}\mathrm{d}\mathrm{a}\mathrm{t}\mathrm{i}\mathrm{o}\mathrm{n}}-{\omega }^2{\mathit{M}}_e\right){\mathit{v}}_e=0.

Equation (32) demonstrates that both thickness variation and form imperfections govern the natural frequency of vibration. So, studying the structure’s buckling, oscillation, and static bending behavior is very interesting because it involves a lot of variables that haven’t been explored in the literature. Furthermore, the computation of the element matrices does not need any shear correction coefficients. This approach guarantees computing efficiency while also improving the accuracy of the mechanical responses of nanoplates.

To use the finite element approach for issue resolution, it is essential to restrict the degrees of freedom at the nodes situated on the border of the plate. If the ith node is situated on the boundary of the plate, where it is subjected to a simply supported condition (represented as S), the constraint condition may be designated as v_{\mathit{\text{z}}\mathrm{b}i}=0,v_{\mathit{\text{z}}\mathrm{s}i}=0. If the node indexed by i is located on the boundary of the plate and is subjected to clamped support (represented by C), the corresponding constraint condition is indicated by v_{\mathit{\text{z}}\mathrm{b}i} =0, v_{\mathit{\text{z}}\mathrm{s}i}= 0, {\left(v_{\mathit{\text{z}}\mathrm{b},x}\right)}_i=0, {\left(v_{\mathit{\text{z}}\mathrm{b},y}\right)}_i=0, {\left(v_{\mathit{\text{z}}\mathrm{s},x}\right)}_i=0, {\left(v_{\mathit{\text{z}}\mathrm{s},y}\right)}_i=0. If the node at position i is located on the free edge (shown as F), there are no limitations or restrictions.

3 Calculation results and discussion

3.1 Verification examples

To demonstrate that the theory put forward in the previous section is valid, this section will show five examples that contrast the results of the nanoplates’ free frequency and critical buckling load with those values obtained through different methods.

Example 1. This instance presents a comparison of the natural frequencies of a (fully simply supported) SSSS plate. The plate has certain geometric and material characteristics. These include the length denoted as a, the thickness indicated as h, where h = a/10 and a/20, respectively. Additionally, the plate exhibits a Young’s modulus of E = 380 GPa, a Poisson’s ratio denoted as V = 0.3, and a mass density represented by \rho = 3960 kg/m³. Two dimensionless elastic foundation parameters are computed as:

(33)\{\begin{array}{c}{K}_1^{}=\frac{k_{\mathrm{w}}{a}^4}{A_0},\\ K_2^{}=\frac{k_{\mathrm{s}}{a}^2}{A_0},\end{array}

where

(34)\begin{array}{rl}& A_0=\frac{E_0{h}^3}{12\left(1-v^2\right)},\\ & E_0=70\mathrm{G}\mathrm{P}\mathrm{a}.\end{array}

The parameter under consideration for calculation and comparison is the dimensionless free frequency of the plate:

(35)K_{\mathrm{f}}={\omega }_{1}h\sqrt{{\rho }_0/E_0},{\rho }_0=2707\mathrm{k}\mathrm{g}/{\mathrm{m}}^3,

Fig.3–Fig.4 show the convergence of the first frequency K_f according to the number of elements; the results in the document [46] are the exact solution. By looking at these results, the following observations are made.

\bulletThe higher the number of elements, the more the calculation results converge and gradually approach the ‘exact’ results in the document [46].

\bulletThe use of a mesh size of 10 × 10, which consists of 100 elements, guarantees the requisite level of precision for the calculations conducted in this study. Consequently, this mesh size will be employed for subsequent cases.

Example 2. The following illustration examines the comparative analysis of the maximum deflection shown by plates when they are subjected to support from an elastic foundation, considering two distinct characteristics. In this case, the plate is subjected to the SSSS boundary condition. The applied stress on the plate is distributed in a sinusoidal fashion, characterized by a magnitude denoted as F₀. The evaluation of plates involves considering their dimensions, which include length (a), width (b), and thickness (h). In this study, we analyze three specific cases. The values of a/h are 5, 10, and 20. The plate’s constituent material has an elastic modulus of 117 GPa and a Poisson’s ratio of 0.33. The equation provided represents the sinusoidal load that is delivered to the plate q=F_0\sin \left({\displaystyle \frac{\text{π}x}{a}}\right)\sin \left({\displaystyle \frac{\text{π}y}{b}}\right) .

The normalization of plate deflection and stress is determined by the following Eq. (36):

(36)\begin{array}{rl}& \overline{w}=\frac{{10}^2{A}_0}{F_0{a}^4}w\left(\frac{a}{2},\frac{b}{2},\frac{h}{2}\right),\\ & {\overline{\sigma }}_x=\frac{1}{100F_0}{\sigma }_x\left(\frac{a}{2},\frac{b}{2},\frac{h}{2}\right),\\ & {\overline{\tau }}_{xy}=\frac{1}{10F_0}{\tau }_{xy}\left(0,0,-\frac{h}{3}\right),\\ & A_0=\frac{E{h}^3}{12\left(1-{\nu }^2\right)}.\end{array}

The elastic foundation is characterized by two parameters that are determined by calculations based on a dimensionless equation:

(37)\{\begin{array}{rl}& K_{\text{w}}^{\ast }=\frac{K_{\text{w}}{a}^4}{A_0},\\ & K_{\text{s}x}^{\ast }=\frac{K_{\text{s}}{a}^2}{A_0},\\ & K_{\text{s}y}^{\ast }=\frac{K_{\text{s}}{b}^2}{A_0}.\end{array}

Tab.1 presents a comparison between the greatest plate deflection observed in this study and the results derived from Navier’s solution [47,48] across various mesh sizes. Tab.1 also demonstrates that increasing the element mesh guarantees convergence of the computation results. Employing a 10-element mesh in every dimension guarantees the necessary degree of precision. Therefore, the ensuing computations use this mesh.

Tab.1 illustrates the comparison between the deflection of a rectangular plate that is supported by an elastic basis. The results obtained by the methodology outlined in this article are determined to be similar to the findings published in earlier studies. This finding provides evidence for the reliability of the computational theory.

Example 3. This example pertains to the calculation of the natural frequency of the structure, denoted as SSSS, with linear variable thickness. The plate has dimensions a = b, and its thickness is given by h_xy(x,y) = h₀(1−{\beta }_x {\displaystyle \frac{x}{a}}). The dimensionaless frequency, which characterizes the plate, is determined by the formula \tilde{\mathrm{\Omega }}=\omega b^2{\displaystyle \frac{\sqrt{\rho h_0/d_0}}{{\text{π}}^2}} (d_0={\displaystyle \frac{E{h}_0^3}{12\left(1-{\nu }^2\right)}}), where E is the Young’s modulus and \rho denotes the mass density. Tab.2 displays the first dimensionless frequency of this plate for various methods, including the current work ({\beta }_k = 0.5), the extended Kantorovich method [49], the spline strip method [50], and the collocation method [51]. It is evident that the data exhibit a high level of agreement.

Example 4. This study focuses on analyzing the critical buckling loads of a plate that is supported by a two-parameter elastic foundation and subjected to an axial compressive stress in the x-direction. The structure has the ratio of its dimensions a/b as 1, and the ratio of its dimensions a/h as 1000. There are two dimensionless elastic foundation coefficients denoted as {\overline{K}}_{\mathrm{w}}={\displaystyle \frac{k_{\mathrm{w}}{a}^4}{D}} and {\overline{K}}_{\mathrm{s}}={\displaystyle \frac{k_{\mathrm{s}}{a}^2}{D}}, where D={\displaystyle \frac{E{h}^3}{12\left(1-{\nu }^2\right)}} and v = 0.3. Tab.3 displays the critical buckling load F_{\mathrm{c}\mathrm{r}}/{\text{π}}^2=F_0{a}^2/\left({\text{π}}^{2}D\right) acquired from the current study, an analytical solution [52], and Green’s functions [53]. The numerical findings indicate a high level of agreement between the data.

Example 5. This last example conducts a buckling analysis of a variable thickness structure. The plate’s thickness is linearly adjusted, and a compressive force is applied in the x-direction. Let us consider a plate characterized by its length, denoted as a, its width, denoted as b, and its thickness, denoted as h. This thickness may be expressed as h₀(1 + {\beta }_x{\displaystyle \frac{x}{a}}), with {\beta }_x= {\displaystyle \frac{h_a}{h_0}}- 1; where h₀ is equal to a divided by 100. At the point x = 0, the thickness of the plate is h₀, whereas at the point x = a, the thickness is h_a. The formula for defining the dimensionless critical buckling load is given as {\tilde{F}}_{\mathrm{c}\mathrm{r}}={\displaystyle \frac{F_0{b}^2}{{\text{π}}^2{D}_e}} (where D_e= {\displaystyle \frac{E{h}_e^3}{12\left(1-{\nu }^2\right)}} and h_e= {\displaystyle \frac{h_0+h_a}{2}}). The data acquired from this particular instance includes the semianalytical solution [54], the Galerkin form [55], and the analog equation method [54,56]. These findings have been compiled in Tab.4. The data exhibits a high degree of similarity among the outcomes.

3.2 Calculation results and discussion

This section gives the numerical calculation results for plates with geometric characteristics where the length of side a is equal to the length of side b. The material parameters of interest in this study are as follows: the elastic stiffness constants c₁₁ and c₃₃ are measured to be 286 and 43.50 GPa, respectively; the piezoelectric coefficient is determined to be 580.3 N/Am; the electromechanical coupling coefficient d₃₃ is found to be 1.57 ×10⁻⁴ N/A²; the piezoelectric charge constant f₁₄ is measured to be 10⁻⁷ N/A; and the mass density of the material is determined to be 5300 kg/m³ [45]. The values assigned to the length scale are {\mathrm{\Gamma }}_0 = {\mathrm{\Gamma }}_1 = {\mathrm{\Gamma }}_2 = 1 nm. A two-parameter elastic foundation with the parameters k_w0 and k_s is what supports the plate.

K_1=\frac{k_{\mathrm{w}0}{a}^4}{c_{11}{h}^3},

K_2=\frac{k_{\mathrm{s}}{a}^2}{c_{11}{h}^3}.

The plate has an initial geometric imperfection as in the theoretical part with the amplitude ratio as shown in Eq. (4) and in Fig.2.

Variation law of plate thickness: h(x) = h₀\left(1+\mathrm{\wp }{\left({\displaystyle \frac{x}{a}}\right)}^n\right)\left(1+\mathrm{\wp }{\left({\displaystyle \frac{y}{b}}\right)}^n\right), where h₀ = a/50 = 10 nm, n shows the change rule of thickness.

In the scenario where n = 0, the plate exhibits a uniform thickness. When n = 1, the plate’s thickness varies linearly. On the other hand, when n = 2, the plate’s thickness demonstrates a nonlinear variation. Additionally, instance \mathrm{\wp } = 0 corresponds to a plate with a constant thickness, denoted as h₀.

Displacement, frequency and critical buckling load of the plate are normalized as follows:

(38)\begin{array}{rl}& N_{\mathrm{d}\mathrm{e}\mathrm{f}\mathrm{l}\mathrm{e}\mathrm{c}\mathrm{t}\mathrm{i}\mathrm{o}\mathrm{n}}=\frac{c_{11}{h}_0}{F_0{a}^2}w,\\ & N_{\mathrm{f}\mathrm{r}\mathrm{e}\mathrm{q}\mathrm{u}\mathrm{e}\mathrm{n}\mathrm{c}\mathrm{y}}={\omega }_i\frac{a^2}{h_0}\sqrt{\frac{\rho }{c_{11}}},\\ & N_{\mathrm{b}\mathrm{u}\mathrm{c}\mathrm{k}}=\frac{N_{cr}{a}^2}{c_{11}{h}_0^3}.\end{array}

To clearly see the influence of the thickness variation factor on the response of the plate, this work provides the following parameters:

(39)\begin{array}{rl}& R_{\omega i}=\frac{{\omega }_i(\mathrm{w}\mathrm{i}\mathrm{t}\mathrm{h}\mathrm{\wp }\ne 0)}{{\omega }_i(\mathrm{w}\mathrm{i}\mathrm{t}\mathrm{h}\mathrm{\wp }=0)},\\ & R_{\mathrm{d}}=\frac{w_{max}(\mathrm{w}\mathrm{i}\mathrm{t}\mathrm{h}\mathrm{\wp }\ne 0)}{w_{max}(\mathrm{w}\mathrm{i}\mathrm{t}\mathrm{h}\mathrm{\wp }=0)},\\ & R_{\mathrm{b}}=\frac{N_{\mathrm{c}\mathrm{r}}(\mathrm{w}\mathrm{i}\mathrm{t}\mathrm{h}\mathrm{\wp }\ne 0)}{N_{\mathrm{c}\mathrm{r}}(\mathrm{w}\mathrm{i}\mathrm{t}\mathrm{h}\mathrm{\wp }=0)}.\end{array}

To clearly see the influence of the variation rule of thickness on the response of the nanoplate, this work gives the parameter n increasing gradually from 0 to 4. The calculation results R_b, R_{\text{ω}1}, R_d are shown in Fig.5–Fig.10. Some comments are as follows.

In the event of a negative value of \mathrm{\wp }, the act of raising n will result in a drop in all of R_b, R_{\text{ω}1}, and R_d. This implies that the plate exhibits the highest level of displacement, the primary natural frequency of oscillation, and a reduced critical buckling load. In the present scenario, the values of R_b, R_{\text{ω}1}, and R_d are all below 1. Additionally, it may be inferred that the plate with varying thickness has a reduced level of stiffness in comparison to the plate with a constant thickness. Furthermore, it should be noted that the value of R_{\text{ω}1} exceeds This phenomenon may be attributed to the variation in plate thickness, which results in a reduction in both the stiffness and mass of the plate. However, it is noteworthy that the oscillation frequency of the plate with variable thickness remains higher than that of the plate with a constant thickness.

In contrast, when \mathrm{\wp } has a positive value, the transformation rule exhibits an inverse relationship compared to the scenario when \mathrm{\wp } possesses a negative value. The plate exhibits enhanced resistance to bending and compression forces in comparison to the scenario where the plate maintains a constant thickness. Due to the fact that R_{\text{ω}1} has a value less than 1, it can be inferred that the frequency of the plate with a variable thickness is greater than the frequency of the plate with a constant thickness. In the scenario when n equals 1, indicating a plate with a linearly changing thickness, the value of R_{\text{ω}1} first lowers to a minimum and thereafter grows. This demonstrates the existence of \mathrm{\wp } for a minimum value of the oscillation frequency of the structure.

Plates exhibiting shape imperfections in various configurations demonstrate consistent transformation rules for R_b and R_d, regardless of the rising values of n. Nevertheless, the degree of variation in parameter R_{\text{ω}1} differs when the shape imperfections exhibit distinct shapes. In particular, in the case of HT, LT1–LT4, the variance of R_{\text{ω}1} progressively rises as n grows. However, the values of the other examples (LT4, LT5, and GT1–GT3) show a declining trend as the variable n increases. It may be inferred from the aforementioned observations that the laws governing thickness variation and form imperfection have a substantial impact on both the stiffness and mass characteristics of a plate that exhibits changing thickness.

The calculated results for R_b, R_{\text{ω}1}, and R_d for plates with varying thickness under various boundary circumstances are shown in Fig.11–Fig.13 and Fig.15. Fig.14 displays the first five natural oscillation frequencies of the plate under two distinct boundary conditions, namely SSSS and SFSF. The presented data illustrates that.

The maximal displacement, the first natural oscillation frequency, and the critical buckling load of the plate vary according to the boundary conditions. In the case of a negative value \mathrm{\wp }, the plate under the SFSF boundary condition has the largest R_b and R_d values, while the plate under the CCCC boundary condition has the smallest R_b and R_d values. In contrast, if \mathrm{\wp } is positive, the opposite occurs.

The transformation rule of variable R_{\text{ω}1}, as dictated by variable \mathrm{\wp }, is contingent upon both the specific nature of the shape imperfection and the prescribed boundary conditions. In the case of a plate exhibiting an imperfect form LT1, subject to the boundary conditions SSSS and SFSF, it is seen that there exists no specific value of \mathrm{\wp } for which the frequency attains its least value. However, for the other boundary conditions, it is possible to determine the minimum frequency value. The frequency value corresponding to the suitable value of \mathrm{\wp } can always be determined for plates labeled as GT1.

The variation of R_b, R_{\text{ω}1}, and R_d depends on the value of {\eta }_{\mathrm{i}\mathrm{m}} shown in Fig.16–Fig.18, it is seen that:

When the value of variable {\eta }_{\mathrm{i}\mathrm{m}} exhibits an upward trend, the transformation rules for variables R_b, R_{\text{ω}1}, and R_d have a high degree of similarity. Specifically, both R_b and R_{\text{ω}1} experience an increase in value over time, with the exception of the GT3 scenario. However, R_d, on the other hand, exhibits a reduction in value. This phenomenon may be elucidated as follows: while the parameter {\eta }_{\mathrm{i}\mathrm{m}} undergoes a rise, the stiffness of the plate also increases. However, the mass of the plate remains constant. Consequently, the displacement of the plate reduces, but the critical load and vibration frequency of the plate experience an increase.

In addition, each type of shape imperfection will have different R_b, R_{\text{ω}1}, and R_d values. This demonstrates that both {\eta }_{\mathrm{i}\mathrm{m}} and shape imperfection types affect the static bending, buckling, and oscillatory bending responses of the plate.

Fig.19–Fig.21 present the results of calculating the critical buckling load, the first natural oscillation frequency, and the maximum displacement of the plate depending on the parameter {\beta }_x of the elastic foundation. The variation of R_b, R_{\text{ω}1}, and R_d depends on the parameter {\beta }_x plotted in Fig.22–Fig.24. These figures lead to the following conclusions:

When the parameter {\beta }_x is raised, the energy of the plate is elevated, resulting in an increase in the plate’s stiffness. As a consequence, the displacement of the plate experiences a reduction, but the frequency and critical load of the plate undergo an increase. The extent to which the frequency, critical buckling load, and maximum displacement vary with increasing {\beta }_x is contingent upon the specific nature of the original form imperfections. The clarity of the degree of change in displacement, frequency, and critical buckling load is less evident in plates with global imperfect type compared to plates with local imperfect type.

The value of R_b consistently exceeds 1, whereas the value of R_d consistently falls below 1, demonstrating that the bending and compression resistance of a plate with varying thickness surpasses that of a plate with a uniform thickness. In the event that the plate exhibits a global imperfect type with two local imperfect kinds (LT4, LT5), the coefficient R_{\text{ω}1}demonstrates a value beyond 1. This finding demonstrates that, in the present scenario, the vibration frequency of the plate exhibiting changing thickness surpasses that of the plate maintaining a constant thickness. In contrast, the reverse phenomenon occurs in the case of HT and LT1–LT3.

4 Conclusions

This paper presents the findings of an investigation into the mechanical behavior of nanoplates with varying thickness, considering the flexomagnetic effect. The study also accounts for imperfections in the initial shape of the nanoplates and incorporates an elastic foundation with coordinate-dependent parameters. The mechanical response of the nanoplates is analyzed in terms of buckling, vibration, and static bending. This study has significant value due to the intricate nature of the calculation procedure and the inherent impossibility of conducting analytical processing. Hence, the present study used the finite element approach. The formulas were derived using an enhanced shear deformation theory. Based on the findings of this study, the following conclusions may be reached.

1) The phenomenon of thickness variation in nanoplates has a significant impact on their mechanical stiffness and overall mass. The stiffness of the plate may either grow or decrease depending on the value of \text{go} when the exponent n is increased.

2) The stiffness of the plate is influenced by both the defect amplitude {\eta }_{\mathrm{i}\mathrm{m}} and the kind of imperfection, whether it is global or local. As the coefficient of stiffness ({\eta }_{\mathrm{i}\mathrm{m}} value) rises, the structural plate exhibits enhanced resistance to compression and bending forces.

3) When the rigidity parameter of the foundation is increased, the plate’s resistance to bending and compressive stresses also improves. However, the extent of the foundation’s influence also depends significantly on the type of imperfection.

The findings of this work contribute to the optimization of the design, fabrication, and use of flexomagnetic nanostructures, particularly in instances where there are imperfections in the initial shape.

References

[1] Van Minh P, Thai L M, Dung N T, Tounsi A, Nhung N T C, Van Thom D. An overview of the flexoelectric phenomenon, potential applications, and proposals for further research directions. International Journal of Mechanics and Materials in Design, 2023, 19(4): 903–925

[2] Liu C, Ke L L, Wang Y S, Yang J, Kitipornchai S. Thermo-electro-mechanical vibration of piezoelectric nanoplates based on the nonlocal theory. Composite Structures, 2013, 106: 167–174

[3] Zhang L L, Liu J X, Fang X Q, Nie G Q. Effects of surface piezoelectricity and nonlocal scale on wave propagation in piezoelectric nanoplates. European Journal of Mechanics. A, Solids, 2014, 46: 22–29

[4] Ke L L, Liu C, Wang Y S. Free vibration of nonlocal piezoelectric nanoplates under various boundary conditions. Physica E, Low-Dimensional Systems and Nanostructures, 2015, 66: 93–106

[5] Li Z, Liu J, Hu B, Wang Y, Shen H. Wave propagation analysis of porous functionally graded piezoelectric nanoplates with a visco-Pasternak foundation. Applied Mathematics and Mechanics, 2023, 44(1): 35–52

[6] Van Minh P, Van Ke T. A Comprehensive study on mechanical responses of non-uniform thickness piezoelectric nanoplates taking into account the flexoelectric effect. Arabian Journal for Science and Engineering, 2023, 48(9): 11457–11482

[7] Zhang L L, Liu J X, Fang X Q, Nie G Q. Size-dependent dispersion characteristics in piezoelectric nanoplates with surface effects. Physica E: Low-Dimensional Systems and Nanostructures, 2014, 57: 169–174

[8] Liu C, Ke L L, Wang Y S, Yang J. Nonlinear vibration of nonlocal piezoelectric nanoplates. International Journal of Structural Stability and Dynamics, 2015, 15(8): 1540013

[9] Kolahchi R, Hosseini H, Esmailpour M. Differential cubature and quadrature-Bolotin methods for dynamic stability of embedded piezoelectric nanoplates based on visco-nonlocal-piezoelasticity theories. Composite Structures, 2016, 157: 174–186

[10] Asemi S R, Farajpour A, Asemi H R, Mohammadi M. Influence of initial stress on the vibration of double-piezoelectric-nanoplate systems with various boundary conditions using DQM. Physica E: Low-Dimensional Systems and Nanostructures, 2014, 63: 169–179

[11] Malekzadeh Fard K, Khajehdehi Kavanroodi M, Malek-Mohammadi H, Pourmoayed A. Buckling and vibration analysis of a double-layer graphene sheet coupled with a piezoelectric nanoplate. Journal of Applied and Computational Mechanics, 2022, 8(1): 129–143

[12] Zhao Z, Zhu J, Chen W. Size-dependent vibrations and waves in piezoelectric nanostructures: A literature review. International Journal of Smart and Nano Materials, 2022, 13(3): 391–431

[13] Du Y M, Lu T, Li X N, Liu Y, Sun W P, Zhang S J, Cheng Z X. High-efficient piezocatalytic hydrogen evolution by centrosymmetric Bi₂Fe₄O₉ nanoplates. Nano Energy, 2022, 104: 107919

[14] Analooei H R, Azhari M, Salehipour H. Thermo-electro-mechanical vibration and buckling analysis of quadrilateral and triangular nanoplates with the nonlocal finite strip method. Mechanics Based Design of Structures and Machines, 2023, 51(3): 1684–1704

[15] Tho N C, Thanh N T, Tho T D, Van Minh P, Hoa L K. Modelling of the flexoelectric effect on rotating nanobeams with geometrical imperfection. Journal of the Brazilian Society of Mechanical Sciences and Engineering, 2021, 43(11): 1–22

[16] LiY DBao RChenW. Buckling of a piezoelectric nanobeam with interfacial imperfection and van der Waals force: Is nonlocal effect really always dominant? Composite Structures, 2018, 194: 357–364

[17] Sahmani S, Aghdam M M. Temperature-dependent nonlocal instability of hybrid FGM exponential shear deformable nanoshells including imperfection sensitivity. International Journal of Mechanical Sciences, 2017, 122: 129–142

[18] Alam M, Mishra S K. Postcritical imperfection sensitivity of functionally graded piezoelectric cylindrical nanoshells using boundary layer solution. Journal of Engineering Mechanics, 2023, 149(9): 04023055

[19] DucD HVan Thom DCongP HVan MinhPNguyenN X. Vibration and static buckling behavior of variable thickness flexoelectric nanoplates. Mechanics Based Design of Structures and Machines, 2022, 51(12): 7102–7130

[20] Wang P, Yuan P, Sahmani S, Safaei B. Surface stress size dependency in nonlinear free oscillations of FGM quasi-3D nanoplates having arbitrary shapes with variable thickness using IGA. Thin-walled Structures, 2021, 166: 108101

[21] Phung V M. Static bending analysis of nanoplates on discontinuous elastic foundation with flexoelectric effect. Journal of Science and Technology, 2022, 17(5): 47–57

[22] Nam V H, Vinh P V, Chinh N V, Thom D V, Hong T T. A new beam model for simulation of the mechanical behaviour of variable thickness functionally graded material beams based on modified first order shear deformation theory. Materials, 2019, 12(3): 404

[23] Van Do T, Hong Doan D, Chi Tho N, Dinh Duc N. Thermal Buckling Analysis of Cracked Functionally Graded Plates. International Journal of Structural Stability and Dynamics, 2022, 22(8): 2250089

[24] Dat P T, Thom D V, Luat D T. Free vibration of functionally graded sandwich plates with stiffeners based on the third-order shear deformation theory. Vietnam Journal of Mechanics, 2016, 38(2): 103–122

[25] Tien D M, Van Thom D, Van Minh P, Tho N C, Doan T N, Mai D N. The application of the nonlocal theory and various shear strain theories for bending and free vibration analysis of organic nanoplates. Mechanics Based Design of Structures and Machines, 2024, 52(1): 588–610

[26] Tuan L T, Dung N T, Van Thom D, Van Minh P, Zenkour A M. Propagation of non-stationary kinematic disturbances from a spherical cavity in the pseudo-elastic cosserat medium. European Physical Journal Plus, 2021, 136(12): 1199

[27] Tran T T, Tran V K, Le P B, Phung V M, Do V T, Nguyen H N. Forced vibration analysis of laminated composite shells reinforced with graphene nanoplatelets using finite element method. Advances in Civil Engineering, 2020, 2020: 1–17

[28] Cong P H, Van Thom D, Duc D H. Phase field model for fracture based on modified couple stress. Engineering Fracture Mechanics, 2022, 269(15): 108534

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

[30] Nanthakumar S S, Lahmer T, Zhuang X, Park H S, Rabczuk T. Topology optimization of piezoelectric nanostructures. Journal of the Mechanics and Physics of Solids, 2016, 94: 316–335

[31] Zhu J, Lv Z, Liu H. Thermo-electro-mechanical vibration analysis of nonlocal piezoelectric nanoplates involving material uncertainties. Composite Structures, 2019, 208: 771–783

[32] Minh P P, Van Do T, Duc D H, Duc N D. The stability of cracked rectangular plate with variable thickness using phase field method. Thin-walled Structures, 2018, 129: 157–165

[33] Nguyen H N, Nguyen T Y, Tran K, Tran T, Nguyen T T, Phan V D, Do T. A finite element model for dynamic analysis of triple-layer composite plates with layers connected by shear connectors subjected to moving load. Materials, 2019, 12(4): 598

[34] Do V T, Pham V V, Nguyen H N. On the development of refined plate theory for static bending behavior of functionally graded plates. Mathematical Problems in Engineering, 2020, 2020: 1–13

[35] Nguyen H N, Hong T T, Vinh P V, Thom D V. An efficient beam element based on Quasi-3D theory for static bending analysis of functionally graded beams. Materials, 2019, 12(13): 2198

[36] Ghasemi H, Park H S, Rabczuk T. A multi-material level set-based topology optimization of flexoelectric composites. Computer Methods in Applied Mechanics and Engineering, 2018, 332: 47–62

[37] Hamdia K M, Ghasemi H, Zhuang X, Alajlan N, Rabczuk T. Sensitivity and uncertainty analysis for flexoelectric nanostructures. Computer Methods in Applied Mechanics and Engineering, 2018, 337: 95–109

[38] Ghasemi H, Park H S, Rabczuk T. A level-set based IGA formulation for topology optimization of flexoelectric materials. Computer Methods in Applied Mechanics and Engineering, 2017, 313: 239–258

[39] Ghasemi H, Park H S, Alajlan N, Rabczuk T. A computational framework for design and optimization of flexoelectric materials. International Journal of Computational Methods, 2020, 17(1): 1850097

[40] Ghasemi H, Park H S, Zhuang X, Rabczuk T. Three-dimensional isogeometric analysis of flexoelectricity with MATLAB implementation. Computers, Materials & Continua, 2020, 65(2): 1157–1179

[41] Soldatos K P. A transverse shear deformation theory for homogeneous monoclinic plates. Acta Mechanica, 1992, 94(3–4): 195–220

[42] Yang W, Liang X, Shen S. Electromechanical responses of piezoelectric nanoplates with flexoelectricity. Acta Mechanica, 2015, 226(9): 3097–3110

[43] Doan D H, Zenkour A M, Van Thom D. Finite element modeling of free vibration of cracked nanoplates with flexoelectric effects. European Physical Journal Plus, 2022, 137(4): 447

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

[45] Zhang N, Zheng S, Chen D. Size-dependent static bending of flexomagnetic nanobeams. Journal of Applied Physics, 2019, 126(22): 223901

[46] Hasani Baferani A, Saidi A R, Ehteshami H. Accurate solution for free vibration analysis of functionally graded thick rectangular plates resting on elastic foundation. Composite Structures, 2011, 93(7): 1842–1853

[47] Sayyad A S, Ghugal Y M. Effects of nonlinear hygrothermomechanical loading on bending of FGM rectangular plates resting on two-parameter elastic foundation using four-unknown plate theory. Journal of Thermal Stresses, 2019, 42(2): 213–232

[48] Mudhaffar I M, Tounsi A, Chikh A, Al-Osta M A, Al-Zahrani M M, Al-Dulaijan S U. Hygro-thermo-mechanical bending behavior of advanced functionally graded ceramic metal plate resting on a viscoelastic foundation. Structures, 2021, 33: 2177–2189

[49] Shufrin I, Eisenberger M. Vibration of shear deformable plates with variable thickness—First-order and higher-order analyses. Journal of Sound and Vibration, 2006, 290(1-2): 465–489

[50] Mlzusawa T. Vibration of rectangular mindlin plates with tapered thickness by the spline strip method. Computers & Structures, 1993, 46(3): 451–463

[51] Takashi M, Jin Y. Application of the collocation method to vibration analysis of rectangular mindlin plates. Computers & Structures, 1984, 18(3): 425–431

[52] Akhavan H, Hashemi S H, Taher H R D, Alibeigloo A, Vahabi S. Exact solutions for rectangular Mindlin plates under in-plane loads resting on Pasternak elastic foundation. Part II: Frequency analysis. Computational Materials Science, 2009, 44(3): 951–961

[53] Lam K Y, Wang C M, He X Q. Canonical exact solutions for Levy-plates on two-parameter foundation using Green’s functions. Engineering Structures, 2000, 22(4): 364–378

[54] Harik I E, Liu X, Ekambaram R. Elastic stability of plates with varying rigidities. Computers & Structures, 1991, 38(2): 161–169

[55] Wittrick W H, Ellen C H. Buckling of tapered rectangular plates in compression. Aeronautical Quarterly, 1962, 13(4): 308–326

[56] Nerantzaki M S, Katsikadelis J T. Buckling of plates with variable thickness—An analog equation solution. Engineering Analysis with Boundary Elements, 1996, 18(2): 149–154

RIGHTS & PERMISSIONS

Higher Education Press