Buckling analysis of nano-scale magneto-electro-elastic plates using the nonlocal elasticity theory

Abstract

Buckling analysis of nonlocal magneto-electro-elastic nano-plate is investigated based on the higher-order shear deformation theory. The in-plane magnetic and electric fields can be ignored for magneto-electro-elastic nano-plates. According to magneto-electric boundary condition and Maxwell equation, the variation of magnetic and electric potentials along the thickness direction of the magneto-electro-elastic plate is determined. To reformulate the elastic theory of magneto-electro-elastic nano-plate, the nonlocal differential constitutive relations of Eringen is applied. Using the variational principle, the governing equations of the nonlocal theory are derived. The relations between local and nonlocal theories are studied by numerical results. Also, the effects of nonlocal parameters, in-plane load directions, and aspect ratio on buckling response are investigated. Numerical results show the effects of the electric and magnetic potentials. These numerical results can be useful in the design and analysis of advanced structures constructed from magneto-electro-elastic materials.

Introduction

In recent years, the study of smart or intelligent materials involving piezoelectric and/or piezomagnetic material structures has attracted many researchers. Smart or intelligent materials have the ability of converting energy among the electric, magnetic, and elastic energies. ¹ In addition, composites made of piezomagnetic/piezoelectric materials show a magnetoelectric effect that is not present in single-phase piezomagnetic or piezoelectric materials. ² Duo to magnetoelectric coupling effects, magneto-electro-elastic (MEE) materials have been widely used in many engineering applications, including the sensor and actuator, robotics, structural health monitoring, vibration control, and medical instruments. For the static deformation of the multilayered MEE plate, the exact closed-form solution is studied by Pan. ³ Some numerical examples in Pan ³ clearly indicate that MEE material possesses special features that may be useful in the design and analysis of smart structures. Also, Pan and Heyliger ⁴ proposed the free vibrations of multilayered MEE plates. Wang et al. ⁵ presented the analytical solution for a three-dimensional (3D) transversely isotropic multilayered MEE simply supported circular plate. Xue et al. ⁶ investigated nonlinear large-deflection model for MEE rectangular plate using the Kirchhoff plate theory. The large deflections of the MEE plate using the meshless local Petrov–Galerkin (MLPG) method is presented by Sladek et al. ⁷ Li ⁸ studied the buckling analysis of MEE plate on elastic foundation. An analytical investigation on buckling and free vibration behavior of Mindlin rectangular MEE nano-plates resting on Pasternak medium via nonlocal elasticity theory has been carried out by Li et al. ⁹ They showed that the normalized frequency of system decreases by increasing the value of electric potential. However, the normalized frequency of the system increases by increasing the value of magnetic potential. Buckling behavior of multilayered MEE plate was investigated by Kiran and Kattimani. ¹⁰ They employed the first-order shear deformation theory and derived the finite element formulation. The analysis of free vibration and biaxial buckling of double-MEE nano-plate-systems subjected to initial external magnetic and electric potentials using nonlocal plate theory was studied by Jamalpoor et al. ¹¹ They supposed that the two nano-plates are bonded with each other using a visco-Pasternak medium. For the analysis of micro- and nano-structures, couple stress theory, modified couple stress theory, and nonlocal elastic theory were used. Buckling analysis of orthotropic protein microtubules under axial and radial compression based on couple stress theory was researched by Beni et al. ¹² Jung et al. ¹³ presented a modified couple stress theory for buckling analysis of sigmoid functionally graded material (S-FGM) nano-plates embedded in Pasternak elastic medium. The nonlocal elastic theory was used in Jung and Han. ¹⁴

From the literature listed above and to the best of authors’ knowledge, however, the buckling analysis under varied aspect ratio of MEE nano-plate on two directions of in-plane load has not been studied. The other buckling studies of MEE plates were confined to x-direction buckling studies. The y-direction buckling studies were limited. Generally, there was no study on the variation of the buckling load in the y-direction due to the change of the aspect ratio and the nonlocal parameters. The in-plane magnetic and electric fields can be ignored for MEE nano-plates. According to magneto-electric boundary condition and Maxwell equation, the variation of magnetic and electric potentials along the thickness direction of the MEE plate is determined. A variational formulation is applied to derive the governing equations of motion and electric, magnetic distribution along the thickness direction of plate by Hamilton principle. In this study, the buckling analysis of a MEE higher-order shear deformable nano-plate is presented in both tabular and graphical forms to investigate the influences of the various parameters. The influences of various aspect ratios, electric and magnetic potentials, and in-plane load directions on the buckling response of MEE nano-plate are investigated. The buckling response obtained by the present method is compared with the results in open literature and a good agreement is acquired.

Higher-order shear deformation theory

The different shape functions $\mathrm{\Psi }(z)$ derived by several researchers are given in Table 1. $\mathrm{\Psi }(z)$ expresses a shape function determining changes in the transverse shear stress and strain distribution along the thickness. In order to obtain accurate results and satisfy the stress-free boundary conditions on the surfaces of the plate, Reddy ¹⁵ proposed a third-order shear deformation plate theory (see Table 1). Touratier ¹⁶ chose transverse strain distribution as a sine function. This theory may be called as sinusoidal shear deformation plate theory. Karama et al. ¹⁷ suggested an exponential variation to investigate the effect of the transverse shear deformation. Mantari et al. ¹⁸ employed the trigonometric shear deformation theory. Hyperbolic shear deformation plate theory was proposed by Mantari and Soares. ¹⁹ In this study, Reddy’s third-order theory with polynomial was selected because it is easy to formulate and expand from the first-order shear deformation theory.

Table 1.

Different shape functions.

Model	Function	$\mathrm{\Psi }(z)$
Reddy 15	Third-order theory with polynomial	$z[1-\frac{4z^2}{3h^2}]$
Touratier 16	Sinusoidal	$\frac{h}{\pi }\sin \left(\frac{\pi z}{h}\right)$
Karama et al. 17	Exponential	$z{e}^{-2{(z/h)}^2}$
Mantari et al. 18	Trigonometric	$\tan \mathrm{mz}-m{\sec }^2\left(\frac{\mathrm{mh}}{2}\right)z,m\ge 0$
Mantari and Soares 19	Hyperbolic	$\sinh \left(\frac{z}{h}\right)e^{\frac{1}{5h}\cosh \left(\frac{z}{h}\right)}+z\left[\frac{\cosh (1/2)+\left[{\sinh }^2(1/2)\right]/5h}{-h}{e}^{\frac{1}{5h}\cosh (1/2)}\right]$

The displacements of a material point located at ( $x,y,z$ ) in the plate may be written as

$$\left\{\begin{array}{c}{u}_x(x,y,z,t)\hfill \\ \multicolumn{1}{c}{u_y(x,y,z,t)}\\ \multicolumn{1}{c}{u_z(x,y,z,t)}\end{array}\right\}=\left\{\begin{array}{c}{u}_x^0(x,y,t)-z\frac{\partial w(x,y,t)}{\partial x}\hfill \\ \multicolumn{1}{c}{u_y^0(x,y,t)-z\frac{\partial w(x,y,t)}{\partial y}}\\ \multicolumn{1}{c}{w(x,y,t)}\end{array}\right\}+\mathrm{\Psi }(z)\left\{\begin{array}{c}{\theta }_x(x,y,t)+\frac{\partial w(x,y,t)}{\partial x}\hfill \\ \multicolumn{1}{c}{{\theta }_y(x,y,t)+\frac{\partial w(x,y,t)}{\partial y}}\\ \multicolumn{1}{c}{0}\end{array}\right\}$$

(1)

in which $u_x^0,u_y^0,w,{\theta }_x,\mathrm{and}{\theta }_y$ are the five unknown displacements of the middle surface, and $\mathrm{\Psi }(z)$ denotes the shape function determining the distribution of the transverse shear stresses and strains along the thickness. By setting $\mathrm{\Psi }(z)=0$ , the displacement field of the classical thin plate theory (CPT) is obtained easily. By setting $\mathrm{\Psi }(z)=z$ , the displacement of the first-order shear deformation plate theory (FSDT) is obtained.

The components of strain tensor, curvature tensor, and rotation vector associated with the displacement field in equation (1) are obtained as

$$\begin{array}{c}\left\{\begin{array}{c}{\epsilon }_{\mathrm{xx}}\hfill \\ \multicolumn{1}{c}{{\epsilon }_{\mathrm{yy}}}\\ \multicolumn{1}{c}{{\gamma }_{\mathrm{xy}}}\end{array}\right\}=\left\{\begin{array}{c}\frac{\partial u_x^0}{\partial x}\hfill \\ \multicolumn{1}{c}{\frac{\partial u_y^0}{\partial y}}\\ \multicolumn{1}{c}{\frac{\partial u_x^0}{\partial y}+\frac{\partial u_y^0}{\partial x}}\end{array}\right\}+z\left\{\begin{array}{c}\frac{\partial {\theta }_x}{\partial x}\hfill \\ \multicolumn{1}{c}{\frac{\partial {\theta }_y}{\partial y}}\\ \multicolumn{1}{c}{\frac{\partial {\theta }_x}{\partial y}+\frac{\partial {\theta }_y}{\partial x}}\end{array}\right\}\hfill \\ \multicolumn{1}{c}{+(\mathrm{\Psi }(z)-z)\left\{\begin{array}{c}\frac{\partial {\theta }_x}{\partial x}+\frac{{\partial }^{2}w}{\partial x^2}\hfill \\ \multicolumn{1}{c}{\frac{\partial {\theta }_y}{\partial y}+\frac{{\partial }^{2}w}{\partial y^2}}\\ \multicolumn{1}{c}{\frac{\partial {\theta }_x}{\partial y}+\frac{\partial {\theta }_y}{\partial x}+2\frac{{\partial }^{2}w}{\partial x\partial y}}\end{array}\right\}}\\ \left\{\begin{array}{c}{\gamma }_{\mathrm{xz}}\hfill \\ \multicolumn{1}{c}{{\gamma }_{\mathrm{yz}}}\end{array}\right\}=\frac{d\mathrm{\Psi }(z)}{\mathrm{dz}}\left\{\begin{array}{c}{\theta }_x+\frac{\partial w}{\partial x}\hfill \\ \multicolumn{1}{c}{{\theta }_y+\frac{\partial w}{\partial y}}\end{array}\right\}\hfill \end{array}$$

(2)

Modeling of the MEE nano-plates

MEE equations

We consider a 3D problem such that all the field variables are functions of the coordinates $(x,y,z)$ , which are coincident with the principal axes of the material symmetry (Figure 1).

Figure 1.

Geometry of the MEE plate.

Denote the symbols $\mathbf{\sigma }$ , $\mathbf{D}$ , and $\mathbf{B}$ as stress, electric displacement, and magnetic induction, respectively. The displacement components along the $x$ , $y$ , and $z$ directions are $u_x^0$ , $u_y^0$ , and $w$ , respectively. The electric potential is denoted as $\varphi$ and the magnetic potential is denoted as $\phi$ . The basic expressions of magneto-electro-elasticity have been given by Parton and Kudryavtsev, ²⁰ Nan, ²¹ and Aboudi. ²² The constitutive equations are

$$\mathbf{\sigma }=\mathbf{C}\mathbf{\gamma }-\mathbf{e}\mathbf{E}-\mathbf{f}\mathbf{H},\mathbf{D}={\mathbf{e}}^{\mathrm{T}}\mathbf{\gamma }+\mathbf{h}\mathbf{E}+\mathbf{g}\mathbf{H},\mathbf{B}={\mathbf{f}}^{\mathrm{T}}\mathbf{\gamma }+\mathbf{g}\mathbf{E}+\mathbf{k}\mathbf{H}$$

(3)

where $\mathbf{\sigma },\mathbf{\gamma },\mathbf{D},\mathbf{E},\mathbf{B}$ , and $\mathbf{H}$ are the vectors of stress, strain, electric displacement, electric field, magnetic induction, and magnetic field, respectively. If the electric vector and magnetic intensity are presented as gradients of the scalar electric and magnetic potentials $\varphi$ and $\phi$ , respectively, Maxwell’s vector equations in the quasi-static approximation are satisfied.

These vectors are defined by

$$\mathbf{\sigma }=\left\{\begin{array}{c}{\sigma }_{\mathrm{xx}}\hfill \\ \multicolumn{1}{c}{{\sigma }_{\mathrm{yy}}}\\ \multicolumn{1}{c}{{\sigma }_{\mathrm{zz}}}\\ \multicolumn{1}{c}{{\sigma }_{\mathrm{xz}}}\\ \multicolumn{1}{c}{{\sigma }_{\mathrm{yx}}}\\ \multicolumn{1}{c}{{\sigma }_{\mathrm{xy}}}\end{array}\right\},\mathbf{\gamma }=\left\{\begin{array}{c}\frac{\partial u_x^0}{\partial x}\hfill \\ \multicolumn{1}{c}{\frac{\partial u_y^0}{\partial y}}\\ \multicolumn{1}{c}{\frac{\partial w}{\partial z}}\\ \multicolumn{1}{c}{\frac{\partial u_x^0}{\partial z}+\frac{\partial w}{\partial x}}\\ \multicolumn{1}{c}{\frac{\partial u_y^0}{\partial z}+\frac{\partial w}{\partial y}}\\ \multicolumn{1}{c}{\frac{\partial u_x^0}{\partial y}+\frac{\partial u_y^0}{\partial x}}\end{array}\right\},\{\mathbf{D},\mathbf{B}\}=\left\{\begin{array}{cc}{D}_x\hfill & B_x\hfill \\ \multicolumn{1}{c}{D_y}& B_y\hfill \\ \multicolumn{1}{c}{D_z}& B_z\hfill \end{array}\right\},\{\mathbf{E},\mathbf{H}\}=\left\{\begin{array}{cc}-\frac{\partial \varphi }{\partial x}\hfill & -\frac{\partial \phi }{\partial x}\hfill \\ \multicolumn{1}{c}{-\frac{\partial \varphi }{\partial y}}& -\frac{\partial \phi }{\partial x}\hfill \\ \multicolumn{1}{c}{-\frac{\partial \varphi }{\partial z}}& -\frac{\partial \phi }{\partial z}\hfill \end{array}\right\}$$

(4)

For a MEE medium whose poling direction coincides with the positive z-axis, the material constant matrices $\mathbf{C},\mathbf{e},\mathbf{f},\mathbf{h},\mathbf{g}$ , and $\mathbf{k}$ have the following forms

$$\begin{array}{c}\mathbf{C}=\left[\begin{array}{cccccc}{C}_{11}\hfill & C_{12}\hfill & C_{13}\hfill & 0\hfill & 0\hfill & 0\hfill \\ \multicolumn{1}{c}{}& C_{22}\hfill & C_{23}\hfill & 0\hfill & 0\hfill & 0\hfill \\ \multicolumn{1}{c}{}& \hfill & C_{33}\hfill & 0\hfill & 0\hfill & 0\hfill \\ \multicolumn{1}{c}{}& \hfill & \hfill & C_{44}\hfill & 0\hfill & 0\hfill \\ \multicolumn{1}{c}{}& \mathrm{sym}\hfill & \hfill & \hfill & C_{55}\hfill & 0\hfill \\ \multicolumn{1}{c}{}& \hfill & \hfill & \hfill & \hfill & C_{66}\hfill \end{array}\right],\mathbf{e}=\left[\begin{array}{ccc}0\hfill & 0\hfill & e_{31}\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & e_{32}\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & e_{33}\hfill \\ \multicolumn{1}{c}{0}& e_{24}\hfill & 0\hfill \\ \multicolumn{1}{c}{e_{15}}& 0\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & 0\hfill \end{array}\right],\mathbf{f}=\left[\begin{array}{ccc}0\hfill & 0\hfill & f_{31}\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & f_{32}\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & f_{33}\hfill \\ \multicolumn{1}{c}{0}& f_{24}\hfill & 0\hfill \\ \multicolumn{1}{c}{f_{15}}& 0\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & 0\hfill \end{array}\right],\hfill \\ \mathbf{h}=\left[\begin{array}{ccc}{h}_{11}\hfill & 0\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& h_{22}\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & h_{33}\hfill \end{array}\right],\mathbf{g}=\left[\begin{array}{ccc}{g}_{11}\hfill & 0\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& g_{22}\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & g_{33}\hfill \end{array}\right],\mathbf{k}=\left[\begin{array}{ccc}{k}_{11}\hfill & 0\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& k_{22}\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & k_{33}\hfill \end{array}\right]\hfill \end{array}$$

(5)

Assuming the absence of the electric charge, body force, and current densities, the equations of motion can be obtained as

$${\sigma }_{\mathrm{ij},j}=\rho \frac{{\partial }^2{u}_i}{\partial t^2},D_{j,j}=0,B_{j,j}=0$$

(6)

with $\rho$ being the density of the material.

Derivation of the governing equations

For the MEE plate, the governing differential equations of motion are derived using Hamilton’s principle which is given as

$$0=\underset{0}{\overset{T}{\int }}(\delta U+\delta (W_1+W_2)-\delta K)\mathrm{dt}$$

(7)

where $\delta U$ is the virtual strain energy, $\delta (W_1+W_2)$ is the virtual work done by applied forces, and $\delta K$ is the virtual kinetic energy.

The strain energy of the MEE plate can be expressed as

$$U=\frac{1}{2}{\displaystyle \underset{\Omega }{\int }{\displaystyle \underset{-h/2}{\overset{h/2}{\int }}\left(\sigma \gamma -DE-BH\right)d\Omega dz}}$$

(8)

Since the MEE plate is thin, the in-plane electric and magnetic field can be deleted, that is, $E_x=E_y=0$ and $H_x=H_y=0$ .

$$U=\frac{1}{2}\underset{\Omega }{\int }\underset{-h/2}{\overset{h/2}{\int }}({\sigma }_{\mathrm{xx}}{\epsilon }_{\mathrm{xx}}+{\sigma }_{\mathrm{yy}}{\epsilon }_{\mathrm{yy}}+{\sigma }_{\mathrm{yz}}{\gamma }_{\mathrm{yz}}+{\sigma }_{\mathrm{zx}}{\gamma }_{\mathrm{zx}}+{\sigma }_{\mathrm{xy}}{\gamma }_{\mathrm{xy}}-D_z{E}_z-B_z{H}_z)d\Omega \mathrm{dz}$$

(9)

Next, we introduce the thickness-integrated stress resultants as

$$M_{\mathrm{ij}}^{(k)}=\underset{-h/2}{\overset{h/2}{\int }}{\sigma }_{\mathrm{ij}}{(z)}^k\mathrm{dz}(k=0,1,2,3)$$

(10)

The kinetic energy of the MEE plate can be described as follows

$$K=\frac{1}{2}\underset{\Omega }{\int }\underset{-h/2}{\overset{h/2}{\int }}[{\left(\frac{\partial u_x}{\partial t}\right)}^2+{\left(\frac{\partial u_y}{\partial t}\right)}^2+{\left(\frac{\partial u_z}{\partial t}\right)}^2]d\Omega$$

(11)

The external work due to transverse loads can be written as

$$W_1=-\underset{\Omega }{\int }{q}_{z}wd\Omega$$

(12)

The potential energy due to in-plane loading is given by

$$W_2=\frac{1}{2}\underset{\Omega }{\int }\left[\begin{array}{cc}\frac{\partial w}{\partial x}\hfill & \frac{\partial w}{\partial y}\hfill \end{array}\right]\left[\begin{array}{cc}{N}_{\mathrm{xi}}\hfill & N_{\mathrm{xyi}}\hfill \\ \multicolumn{1}{c}{N_{\mathrm{xyi}}}& N_{\mathrm{yi}}\hfill \end{array}\right]\left[\begin{array}{c}\frac{\partial w}{\partial x}\hfill \\ \multicolumn{1}{c}{\frac{\partial w}{\partial y}}\end{array}\right]d\Omega (i=m,e,q)$$

(13)

$N_x$ , $N_{\mathrm{xe}}$ , $N_{\mathrm{xm}}$ , $N_y$ , $N_{\mathrm{ye}}$ , and $N_{\mathrm{ym}}$ are, respectively, the mechanical, electric, and magnetic forces for the $x-$ and $y-$ directions with

$$N_x=P,N_{\mathrm{xe}}=e_{31}{\mathrm{V}}_0,N_{\mathrm{xm}}=f_{31}{\Omega }_0$$

(14)

$$N_y=\lambda P,N_{\mathrm{ye}}=e_{31}{\mathrm{V}}_0,N_{\mathrm{ym}}=f_{31}{\Omega }_0$$

(15)

in which $P$ is the mechanical load along $x-$ direction and $\lambda$ is the lateral load parameter.

Nonlocal elastic theory

The stress at a point depends only on the strain at that point in classical local elasticity theories. While in nonlocal elasticity theories, it is assumed that the stress at a point depends on the strains at all the points of the continuum. In other words, according to this nonlocal theory, strain at a point depends on both stress and spatial derivatives of the stress at that point. The nonlocal constitutive behavior of a Hookean solid (Eringen^23,24) is expressed by the following differential constitutive relation

$$(1-\mu {\nabla }^2){\sigma }_{\mathrm{ij}}=t_{\mathrm{ij}}$$

(16)

where $t_{\mathrm{ij}}$ is the local stress tensor, ${\sigma }_{\mathrm{ij}}$ in the nonlocal stress tensor, and the nonlocal parameter $\mu$ is defined by

$$\mu =e_0^2{\overline{a}}^2$$

(17)

in which $e_0$ is the material constant which is defined by the experiment and $\overline{a}$ is the internal characteristic length. ²⁵

The relations between stress resultants in local and nonlocal theory are determined by integrating equation (16) through the plate thickness

$$\mathcal{L}\left(M_{ij}^{(0)}\right)=M_{ij}^{(0)L},\mathcal{L}\left(M_{ij}^{(1)}\right)=M_{ij}^{(1)L},\mathcal{L}\left(M_{ij}^{(2)}\right)=M_{ij}^{(2)L}$$

(18)

where $\mathcal{L}=1-\mu {\nabla }^2$ and

$$\left\{\begin{array}{c}{M}_{\alpha \beta }^{(0)},M_{\alpha \beta }^{(0)L}\hfill \\ \multicolumn{1}{c}{M_{\alpha \beta }^{(1)},M_{\alpha \beta }^{(1)L}}\\ \multicolumn{1}{c}{M_{\alpha \beta }^{(2)},M_{\alpha \beta }^{(2)L}}\end{array}\right\}=\underset{-h/2}{\overset{h/2}{\int }}\{{\sigma }_{\alpha \beta },t_{\alpha \beta }\}\left\{\begin{array}{c}1\hfill \\ \multicolumn{1}{c}{z}\\ \multicolumn{1}{c}{z^3}\end{array}\right\}\mathrm{dz}$$

(19)

$$\left\{\begin{array}{c}{M}_{\alpha z}^{(0)},M_{\alpha z}^{(0)L}\hfill \\ \multicolumn{1}{c}{M_{\alpha z}^{(2)},M_{\alpha z}^{(2)L}}\end{array}\right\}=\underset{-h/2}{\overset{h/2}{\int }}\{{\sigma }_{\alpha z},t_{\alpha z}\}\left\{\begin{array}{c}1\hfill \\ \multicolumn{1}{c}{z^2}\end{array}\right\}\mathrm{dz}$$

(20)

where $\alpha$ and $\beta$ take the symbols $x$ , $y$ . The superscript L denotes the symbol in local third-order shear deformation theory and $h$ is the thickness of the plate.

In general, differential operator ∇ in equation (16) is the 3D Laplace operator. For 2D problems, the operator ∇ may be reduced to 2D one. Thus, the linear differential operator $\mathcal{L}$ becomes

$$\overline{\mathcal{L}}=1-\mu (\frac{{\partial }^2}{\partial x^2}+\frac{{\partial }^2}{\partial y^2})$$

(21)

It is clear that the operator $\overline{\mathcal{L}}$ is independent of $z$ direction.

Finally, using Hamilton’s principles (equation (7)) lead to the following governing equations

$$\delta u_x^0:M_{\mathrm{xx},x}^{(0)}+M_{\mathrm{xy},y}^{(0)}-\overline{\mathcal{L}}\left(I_0{\stackrel{··}{u}}_x^0\right)=0$$

(22a)

$$\delta u_y^0:M_{\mathrm{xy},x}^{(0)}+M_{\mathrm{yy},y}^{(0)}-\overline{\mathcal{L}}\left(I_0{\stackrel{··}{u}}_y^0\right)=0$$

(22b)

$$\begin{array}{c}\delta w:M_{\mathrm{zx},x}^{(0)}+M_{\mathrm{zy},y}^{(0)}-c_2(M_{\mathrm{zx},x}^{(2)}+M_{\mathrm{zy},y}^{(2)})\hfill \\ \multicolumn{1}{c}{+c_1(M_{\mathrm{xx},\mathrm{xx}}^{(2)}+M_{\mathrm{yy},\mathrm{yy}}^{(2)}+2M_{\mathrm{xy},\mathrm{xy}}^{(2)})}\\ -\overline{\mathcal{L}}[-q_z-N_0+I_0\stackrel{··}{w}+c_1{I}_4({\stackrel{··}{\theta }}_{x,x}+{\stackrel{··}{\theta }}_{y,y})-c_1^2{I}_6({\stackrel{··}{\theta }}_{x,x}+{\stackrel{··}{\theta }}_{y,y}+{\stackrel{··}{w}}_{,\mathrm{xx}}+{\stackrel{··}{w}}_{,\mathrm{yy}})]=0\hfill \end{array}$$

(22c)

$$\delta {\theta }_x:M_{xx,x}^{(1)}+M_{xy,y}^{(1)}-c_1\left(M_{xx,x}^{(2)}+M_{xy,y}^{(2)}\right)-M_{zx}^{(0)}+c_2{M}_{zx}^{(2)}-\overline{\mathcal{L}}\left[I_2{\stackrel{··}{\theta }}_x-c_1{I}_4\left({\stackrel{··}{\theta }}_x+{\stackrel{··}{w}}_{,x}\right)-c_1\left(I_4{\stackrel{··}{\theta }}_x-c_1{I}_6\left({\stackrel{··}{\theta }}_x+{\stackrel{··}{w}}_{,x}\right)\right)\right]=0$$

(22d)

$$\delta {\theta }_y:M_{xy,x}^{(1)}+M_{yy,y}^{(1)}-c_1\left(M_{xy,x}^{(2)}+M_{yy,y}^{(2)}\right)-M_{zy}^{(0)}+c_2{M}_{zy}^{(2)}-\overline{\mathcal{L}}\left[I_2{\stackrel{··}{\theta }}_y-c_1{I}_4\left({\stackrel{··}{\theta }}_y+{\stackrel{··}{w}}_{,y}\right)-c_1\left(I_4{\stackrel{··}{\theta }}_y-c_1{I}_6\left({\stackrel{··}{\theta }}_y+{\stackrel{··}{w}}_{,y}\right)\right)\right]=0$$

(22e)

$$\frac{\partial D_z}{\partial z}=0$$

(22f)

$$\frac{\partial B_z}{\partial z}=0$$

(22g)

where

$$\begin{array}{c}{N}_0=(N_{\mathrm{xm}}+N_{\mathrm{xe}}+N_{\mathrm{xq}})\frac{{\partial }^{2}w}{\partial x^2}+(N_{\mathrm{ym}}+N_{\mathrm{ye}}+N_{\mathrm{yq}})\frac{{\partial }^{2}w}{\partial y^2}\hfill \\ I_k=\underset{-h/2}{\overset{h/2}{\int }}\rho {\left(z\right)}^k\mathrm{dz}(k=0,2,4,6)\hfill \\ c_2=3c_1\hfill \\ c_1=\frac{4}{3h^2}\hfill \end{array}$$

in which $\rho$ is the mass density. From these, by setting $\overline{\mathcal{L}}=1$ , the equations of motion of the local third-order shear deformation theory are obtained.

Substituting equation (4) into equations (22f) and (22g), the two algebraic equations can be written as

$$e_{31}\left(\frac{\partial {\theta }_x}{\partial x}-c_2\left(\frac{\partial {\theta }_x}{\partial x}+\frac{{\partial }^{2}w}{\partial x^2}\right)\right)+e_{31}\left(\frac{\partial {\theta }_y}{\partial y}-c_2\left(\frac{\partial {\theta }_y}{\partial y}+\frac{{\partial }^{2}w}{\partial y^2}\right)\right)-h_{33}\frac{{\partial }^2\varphi }{\partial z^2}-g_{33}\frac{{\partial }^2\phi }{\partial z^2}=0$$

(23)

$$f_{31}\left(\frac{\partial {\theta }_x}{\partial x}-c_2\left(\frac{\partial {\theta }_x}{\partial x}+\frac{{\partial }^{2}w}{\partial x^2}\right)\right)+f_{31}\left(\frac{\partial {\theta }_y}{\partial y}-c_2\left(\frac{\partial {\theta }_y}{\partial y}+\frac{{\partial }^{2}w}{\partial y^2}\right)\right)-g_{33}\frac{{\partial }^2\varphi }{\partial z^2}-k_{33}\frac{{\partial }^2\phi }{\partial z^2}=0$$

(24)

By adopting Crammer’s rule, one can be obtained as

$$\frac{{\partial }^2\varphi }{\partial z^2}=M_1\Delta ,\frac{{\partial }^2\phi }{\partial z^2}=M_2\Delta$$

(25)

where

$$\begin{array}{c}{M}_1=\frac{k_{33}{e}_{31}-g_{33}{f}_{31}}{h_{33}{k}_{33}-g_{33}^2},M_2=\frac{h_{33}{f}_{31}-g_{33}{e}_{31}}{h_{33}{k}_{33}-g_{33}^2}\hfill \\ \Delta =\frac{\partial {\theta }_x}{\partial x}+\frac{\partial {\theta }_y}{\partial y}-c_2(\frac{\partial {\theta }_x}{\partial x}+\frac{\partial {\theta }_y}{\partial y}+\frac{{\partial }^{2}w}{\partial x^2}+\frac{{\partial }^{2}w}{\partial y^2})\hfill \end{array}$$

(26)

It can be derived from equations (25) and (26) that

$$\frac{\partial \varphi }{\partial z}=M_{1}z\Delta +{\varphi }_0(x,y),\frac{\partial \phi }{\partial z}=M_{2}z\Delta +{\phi }_0(x,y)$$

(27)

and

$$\varphi =\frac{M_1\Delta }{2}(z^2-\frac{h^2}{4})+\frac{V_0}{h}z+\frac{V_0}{2},\phi =\frac{M_2\Delta }{2}(z^2-\frac{h^2}{4})+\frac{{\Omega }_0}{h}z+\frac{{\Omega }_0}{2}$$

(28)

where the electric and magnetic boundary conditions are determined as $\varphi (h/2)=V_0$ , $\varphi (-h/2)=0$ and $\phi (h/2)={\Omega }_0$ , $\phi (-h/2)=0$ .

Buckling results of nonlocal third-order shear deformation theory

Here, analytical solutions for buckling response of simply supported MEE nano-plates are investigated using the nonlocal third-order shear deformation theory to illustrate the magnetoelectric effects on buckling loads. In the case of simply supported boundary conditions, the Navier solution can be easily obtained. According to the Navier method, the displacements at the middle surface of the plane ( $z=0$ ) are represented in double-Fourier series as

$$\left\{\begin{array}{c}\begin{array}{c}{u}_x^0(x,y,t)\hfill \\ \multicolumn{1}{c}{u_y^0(x,y,t)}\\ \multicolumn{1}{c}{w(x,y,t)}\\ \multicolumn{1}{c}{{\theta }_x(x,y,t)}\end{array}\hfill \\ \multicolumn{1}{c}{{\theta }_y(x,y,t)}\end{array}\right\}=\sum _{m,n=1}^{\infty }\left\{\begin{array}{c}\begin{array}{c}{U}_{\mathrm{mn}}(t){\Lambda }_1\hfill \\ \multicolumn{1}{c}{V_{\mathrm{mn}}(t){\Lambda }_2}\\ \multicolumn{1}{c}{W_{\mathrm{mn}}(t){\Lambda }_3}\\ \multicolumn{1}{c}{X_{\mathrm{mn}}(t){\Lambda }_1}\end{array}\hfill \\ \multicolumn{1}{c}{Y_{\mathrm{mn}}(t){\Lambda }_2}\end{array}\right\}{e}^{i\omega t}$$

(29)

where ${\Lambda }_1=\cos \xi x\sin \eta y$ , ${\Lambda }_2=\sin \xi x\cos \eta y$ , and ${\Lambda }_3=\sin \xi x\sin \eta y$ , in which, $\xi =m\pi /a$ and $\eta =n\pi /b$ .

For the simply supported MEE nano-plate, we have the following boundary conditions

$$u_y^0=w={\theta }_y=M_{\mathrm{xx}}^{(0)}=M_{\mathrm{xx}}^{(1)}=M_{\mathrm{xx}}^{(2)}=0,\mathrm{at}x=0,a$$

(30)

$$u_x^0=w={\theta }_x=M_{\mathrm{yy}}^{(0)}=M_{\mathrm{yy}}^{(1)}=M_{\mathrm{yy}}^{(2)}=0,\mathrm{at}y=0,b$$

(31)

By substituting equation (29) into equations (22a)–(22e), we obtain the following matrix form

$$\left[\mathbf{K}\right]\left\{\Delta \right\}+\left[\mathbf{M}\right]\left\{\stackrel{··}{\Delta }\right\}=\left\{Q\right\}$$

(32)

where $\{\Delta \}=\{U_{\mathrm{mn}},V_{\mathrm{mn}},W_{\mathrm{mn}},X_{\mathrm{mn}},Y_{\mathrm{mn}}\}$ , $[\mathbf{K}]$ is the stiffness matrix (see Appendix 1), $[\mathbf{M}]$ is the mass matrix, and $\{Q\}$ is the transverse load vector. The transverse load $q_z$ is also represented in the double-Fourier sine series as

$$q_z(x,y)=\sum _{m=1}^{\infty }\sum _{n=1}^{\infty }{Q}_{\mathrm{mn}}\sin \xi x\sin \eta y$$

(33)

For static bending analysis, equation (32) takes $[\mathbf{M}]=0$ and for free vibration analysis, set $\{Q\}$ to zero in equation (32).

For buckling analysis, the matrix form is as follows

$$\left[\mathbf{K}\right]\left\{\Delta \right\}-{\lambda }_{\mathrm{mn}}\left[\mathbf{G}\right]\left\{\Delta \right\}=(\left[\mathbf{K}\right]-{\lambda }_{\mathrm{mn}}\left[\mathbf{G}\right])\left\{\Delta \right\}=\left\{0\right\}$$

(34)

where ${\lambda }_{\mathrm{mn}}$ are the eigenvalues and $[\mathbf{G}]=0$ , except ${\mathbf{G}}_{33}={\xi }^2+(N_y/N_x){\eta }^2$ . Equation (34) defines a standard eigenvalue problem, with ${\lambda }_{\mathrm{mn}}=N_{\mathrm{cr}}(m,n)$ . The eigenvalues ${\lambda }_{\mathrm{mn}}$ are calculated by setting the determinant of the coefficient matrix in equation (34) to zero. The smallest eigenvalue is the critical buckling load (Figure 2).

Figure 2.

Loading conditions of FGM nano-plate for uniaxial compression (a) $N_x$ and (b) $N_y$ .

Numerical results

In order to confirm the present higher-order plate theory, some convergence analyses are performed to verify the accuracy of the derived solutions and formulations against the references.

Validation

First, the buckling loads for single MEE plate are calculated and compared with the results in Li et al. ⁹ As shown in Table 2, the present buckling loads are in good agreement with those studied in Li et al. ⁹ Due to the effect of the higher-order shear deformation, in the case of side-to-thickness ratio ( $a/h=10$ ), the difference between the present results and Li et al. ⁹ is found.

Table 2.

Buckling loads for MEE plate ( $V_0={\Omega }_0=0,m=n=1$ ).

$P_{\mathrm{cr}}$	Side-to-thickness ratio ( $a/h$ )
	10		100		1000
	Present	Liv et al. 9	Present	Li et al. 9	Present	Li et al. 9
$\lambda =0$	2.978643	2.9747	3.27604	3.2760	3.279316	3.2794
$\lambda =0.5$	1.985762	1.9831	2.18403	2.1840	2.186211	2.1862
$\lambda =-0.5$	5.957285	5,9494	6.55208	6.5520	6.558632	6.5587

In second analysis, as can be demonstrated in Figure 3, for MEE nano-plate, the dimensionless buckling loads calculated for various values of the nonlocal parameter from this study are compared with the results in Li et al. ⁹ The results of this study showed good agreement.

Figure 3.

Variation of buckling load, $P_{\mathrm{cr}}$ , versus nonlocal parameter, $e_{0}a$ , for MEE nano-scale plate under different lateral load parameter $\lambda$ ( $V_0={\Omega }_0=0,m=n=1)$ .

Third, the effects of electric and magnetic potential on the buckling load for MEE plate are illustrated in Figure 4. It is shown that the buckling load decreases linearly with an increase in the value of electric potential. Contrary to the type of electric potential shown in Figure 4, the buckling load increases with the increase in magnetic potential.

Figure 4.

Variation of buckling load, $P_{\mathrm{cr}}$ , versus electrical and magnetic potential, $V_0,{\Omega }_0$ , for MEE plate under different side-to-thickness ratio $a/h$ ( $m=n=1)$ .

Parameter study

Parameter studies of buckling analysis for MEE nano-plate are presented. The MEE nano-plates made of BaTiO₃ (piezoelectric material) as the inclusions and CoFe₂O₄ (piezomagnetic material) as the matrix are considered. The material properties of the MEE nano-plate are as follows: ⁹

$$\begin{array}{c}{C}_{11}=22.6\times {10}^{10}\mathrm{N}/{\mathrm{m}}^2,C_{12}=12.4\times {10}^{10}\mathrm{N}/{\mathrm{m}}^2,\hfill \\ \multicolumn{1}{c}{C_{22}=21.6\times {10}^{10}\mathrm{N}/{\mathrm{m}}^2,C_{44}=4.4\times {10}^{10}\mathrm{N}/{\mathrm{m}}^2,}\\ e_{31}=-2.2\mathrm{C}/{\mathrm{m}}^2,f_{31}=290.2\mathrm{N}/(\mathrm{Am}),\hfill \\ h_{33}=63.5\times {10}^{-10}{\mathrm{C}}^2/(\mathrm{N}{\mathrm{m}}^2),\hfill \\ \multicolumn{1}{c}{g_{33}=27.375\times {10}^{-10}\mathrm{Ns}/(\mathrm{VC}),}\\ k_{33}=835,000\times {10}^{-10}\mathrm{N}{\mathrm{s}}^2/{\mathrm{C}}^2\hfill \end{array}$$

(35)

In the research of MEE nano-plate, actually, the material properties must be size-dependent and should be determined by molecular dynamic simulation or experiments. However, currently, such studies on the MEE nano-materials are lacking. Therefore, in order to understand the magneto-electro-mechanical responses of the MEE nano-plate, we choose the material properties of the macroscopic MEE materials for the case study. The length and the width of the nano-plate are a = b = 10 nm. The thickness of the nano-plate is taken as h = 1 nm. The results are presented in the nondimensional form using equation (36)

$$P_{\mathrm{cr}}=\frac{P{a}^2}{C_{11}{h}^3}$$

(36)

The effects of aspect ratio (b/a) on nondimensional buckling load of MEE nano-scale plate under uniaxial compression $N_x$ and $N_y$ are shown in Figures 5–9. The side-to-thickness ratio ( $a/h$ ) is assumed to be 10 and mode number assumed to be 1 except for Figure 8. In addition, the electrical potential $V_0=50$ and magnetic potential ${\Omega }_0=5$ are considered. In the case of nonlocal theory, as shown in Figure 8, the graph is overlapped due to the change of buckling mode as the respect ratio b/a < 0.4.

Figure 5.

Variation of buckling load, $P_{\mathrm{cr}}$ , versus aspect ratio, (b/a), for MEE nano-scale plate under different nonlocal parameter $\mu$ ( $V_0={\Omega }_0=0,m=n=1,a/h=10)$ .

Figure 6.

Variation of buckling load, $P_{\mathrm{cr}}$ , versus aspect ratio, (b/a), for MEE nano-scale plate under different nonlocal parameter $\mu$ ( $V_0=0,{\Omega }_0=5,m=n=1,a/h=10)$ .

Figure 7.

Variation of buckling load, $P_{\mathrm{cr}}$ , versus aspect ratio, (b/a), for MEE nano-scale plate under different nonlocal parameter $\mu$ ( $V_0=50,{\Omega }_0=0,m=n=1,a/h=10)$ .

Figure 8.

Variation of buckling load, $P_{\mathrm{cr}}$ , versus aspect ratio, (b/a), for MEE nano-scale plate under different mode number $m,n$ ( $V_0=50,{\Omega }_0=0,\mu =1.00,a/h=10)$ .

Figure 9.

Variation of buckling load, $P_{\mathrm{cr}}$ , versus aspect ratio, (b/a), for MEE nano-scale plate under different nonlocal parameter $\mu$ ( $V_0=50,{\Omega }_0=5,m=n=1,a/h=10)$ .

To account for the effect of nonlocal parameter on responses of MEE nano-scale plate, Figure 5 plots the nondimensional buckling load with respect to the nonlocal parameter for a simply supported MEE nano-scale plate with $m=n=1$ , $V_0={\Omega }_0=0$ , and $a/h=10$ . Figure 5 clearly shows the diminishing effect of aspect ratio (b/a) on buckling loads, the effect being negligible for aspect ratio larger than 2. The increasing value of the nonlocal parameter leads to a decrease in the magnitude of buckling load. The increasing value of aspect ratio leads to a decrease in the buckling load. As expected, the effect of nonlocal parameter decreases the buckling load.

The buckling load for the magnetic potential ${\Omega }_0=5$ is shown in Figure 6. Magnetic potential has the effect of increasing the buckling load. Overall buckling load was increased by 5 times. It was found that the magnetic potential plays a role in increasing the stiffness of the MEE nano-scale plate. As expected, increasing value of the aspect ratio and the nonlocal parameter leads to a decrease in the magnitude of buckling load. The increasing value of the nonlocal parameter plays a role in decreasing the stiffness of the MEE nano-scale plate.

In Figure 7, buckling loads for the electrical potential $V_0=50$ are investigated. As the nonlocal parameter increases, the buckling load changes at aspect ratio (b/a) = 0.2. It is considered that the mixed effect of the electrical potential and the nonlocal parameter caused the special situation. Especially, when the aspect ratio (b/a) <0.2, the sign of the buckling load changed. It is considered that the mixed effect of the nonlocal parameter and the effect of electrical potential changed the sign of stiffness. These results are considered to have no physical meaning and are considered to be a weak point of numerical interpretation. In the future, it will be necessary to clarify the exact physical limitations accompanied with experimental studies. Until experimental identification, it is considered that only aspect ratio (b/a)> 0.2 should be considered.

In the case of electrical potential $V_0=50$ , the buckling load according to the mode number is presented in Figure 8. As shown in Figure 7, buckling load changes when the aspect ratio is less than 0.3. As expected, if the nonlocal parameter is not taken into account, the buckling load overlaps when the aspect ratio is less than 0.8. In addition, buckling loads exceeding the physical limit due to the influence of electrical potential were shown when nonlocal parameters were considered.

The effects of aspect ratio (b/a) on nondimensional buckling load of MEE nano-scale plate under uniaxial compression $N_x$ are shown in Figure 9. The side-to-thickness ratio ( $a/h$ ) is assumed to be 10 and mode number assumed to be 1 except for Figure 9. In addition, the electrical potential and magnetic potential, $V_0=50,{\Omega }_0=5$ , are considered. The effect of increasing the stiffness of the MEE nano-scale plate by the magnetic potential is overwhelmingly greater than the effect of reducing the stiffness by the electrical potential. Therefore, even if the aspect ratio is less than 1, the parabolic change of the buckling load does not occur.

In Figure 10, the case of the MEE nano-scale plate subjected to y-directional compression, the parabolic change of the buckling load occurred when the aspect ratio was less than 1, even when the magnetic potential and electrical potential were zero. Numerical (physical) limitations of the MEE nano-scale plate could be presented according to the variation of nonlocal parameter.

Figure 10.

Variation of buckling load, $P_{\mathrm{cr}}$ , versus aspect ratio, (b/a), for MEE nano-scale plate under different nonlocal parameter $\mu$ ( $V_0={\Omega }_0=0,m=n=1,a/h=10)$ .

Figure 11 shows the buckling load according to the change of buckling mode. As can be expected, when the aspect ratio is less than 1, the buckling load changes in parabolic shape are shown. When the nonlocal parameter was considered, no parabolic change was observed in the mode (m = 2, n = 1).

Figure 11.

Variation of buckling load, $P_{\mathrm{cr}}$ , versus aspect ratio, (b/a), for MEE nano-scale plate under different mode number $m,n$ ( $V_0={\Omega }_0=0,\mu =0,1.0a/h=10)$ .

Figures 12–14 show the buckling load in consideration of magnetic potential and electrical potential. In Figure 12, the sign of the buckling load changes, which is analyzed as the effect of decreasing the stiffness by the electrical potential, as in the case of $N_x$ . Figure 12 clearly shows that the increasing value of the nonlocal parameter leads to a decrease in the magnitude of buckling load. The effect of aspect ratio is decreased when the aspect ratio is 1, and then increased again. As the nonlocal parameter increases, the buckling load changes at aspect ratio (b/a) = 0.25. It is considered that the mixed effect of the electrical potential and the nonlocal parameter caused the special situation.

Figure 12.

Variation of buckling load, $P_{\mathrm{cr}}$ , versus aspect ratio, (b/a), for MEE nano-scale plate under different nonlocal parameter $\mu$ ( $V_0=50,{\Omega }_0=0,m=n=1,a/h=10)$ .

Figure 13.

Variation of buckling load, $P_{\mathrm{cr}}$ , versus aspect ratio, (b/a), for MEE nano-scale plate under different nonlocal parameter $\mu$ ( $V_0=0,{\Omega }_0=5,m=n=1,a/h=10)$ .

Figure 14.

Variation of buckling load, $P_{\mathrm{cr}}$ , versus aspect ratio, (b/a), for MEE nano-scale plate under different nonlocal parameter $\mu$ ( $V_0=50,{\Omega }_0=5,m=n=1,a/h=10)$ .

For the magnetic potential ${\Omega }_0=5$ , the effects of aspect ratio (b/a) on nondimensional buckling load of MEE nano-scale plate under uniaxial compression $N_y$ are shown in Figure 13. Magnetic potential has the effect of increasing the buckling load. Overall buckling load was increased by 5 times. It was found that the magnetic potential plays a role in increasing the stiffness of the MEE nano-scale plate.

In Figure 14, the electrical potential and magnetic potential, $V_0=50,{\Omega }_0=5$ , are considered. The effect of increasing the stiffness of the MEE nano-scale plate by the magnetic potential is overwhelmingly greater than the effect of reducing the stiffness by the electrical potential. Unlike the case of $N_x$ , even if the aspect ratio is less than 0.25, the parabolic change of the buckling load occurred.

Figure 15 illustrates the buckling load according to the buckling mode. As in the case of $N_x$ , parabolic shape buckling load change is shown when the aspect ratio is smaller than 1. When the nonlocal parameter was considered, no parabolic change occurred in the mode (m = 2, n = 1).

Figure 15.

Variation of buckling load, $P_{\mathrm{cr}}$ , versus aspect ratio, (b/a), for MEE nano-scale plate under different mode number $m,n$ ( $V_0=50,{\Omega }_0=5,\mu =0,1.0,a/h=10)$ .

Conclusion

In this article, the buckling analysis of MEE nano-plate is investigated based on nonlocal theory. The in-plane magnetic and electric fields can be ignored for MEE nano-plates. According to magneto-electric boundary condition and Maxwell equation, the variation of magnetic and electric potentials along the thickness direction of the MEE plate is determined. From the numerical results, some conclusions can be drawn:

1. The buckling load decreases with increasing nonlocal parameter for a MEE nano-plate.

2. For a MEE nano-plate, the buckling load decreases linearly with electric potential, and increases with magnetic potential.

3. The effect of varying the aspect ratio according to the direction of the in-plane compressive load is different. Therefore, it is considered that the change of the direction and the aspect ratio of the in-plane compressive load should be considered in two directions, respectively.

4. The comparison of the buckling load between the case with and without the magnetic potential and the electrical potential shows the difference of the buckling load according to the buckling mode in the case of the nonlocal MEE nano-plate under N_y.

The results given in this study can be used as benchmark for buckling behavior of other nano-scale MEE plate. Furthermore, the presented theory may be extended to other cases of MEE structures such as MEE shells on elastic foundation.