Influence of flexoelectric,small-scale,surface and residual stress on the nonlinear vibration of sigmoid,exponential and power-law FG Timoshenko nano-beams

Abstract

This research deals with the nonlinear vibration of the functionally graded nano-beams based on the nonlocal elasticity theory considering surface and flexoelectric effects. The flexoelectric functionally graded nano-beam is resting on nonlinear Pasternak foundation. Cubic nonlinearity is assumed for foundation. It is assumed that the material properties of the nano-beam change continuously along the thickness direction according to different patterns of material distribution. In order to include coupling of strain gradients and electrical polarizations in equation of motion, the nonlocal, nonclassical nano-beam model containing flexoelectric effect is employed. In addition, the effects of surface elasticity, di-electricity, and piezoelectricity as well as bulk flexoelectricity are accounted in constitutive relations. The governing equations of motion are derived using Hamilton principle based on first shear deformation beam theory and the nonlocal strain gradient elasticity theory considering residual surface stresses. The differential quadrature method is used to calculate nonlinear natural frequency of flexoelectric functionally graded nano-beam as well as nonlinear vibrational mode shape. After validation of the present numerical results with those results available in literature, full numerical results are presented to investigate the influence of important parameters such as flexoelectric coefficients of the surface and bulk, residual surface stresses, nonlocal parameter, length scale effects (strain gradient parameter), cubic nonlinear Winkler and shear coefficients, power gradient index of functionally graded material, and geometric dimensions on the nonlinear vibration behaviors of flexoelectric functionally graded nano-beam. The numerical results indicate that, considering the flexoelectricity leads to the decrease of the bending stiffness of the flexoelectric functionally graded nano-beams.

Keywords

Flexoelectric surface effects residual surface stresses functionally graded beam nonlinear vibration sigmoid and power

Introduction

The beams are one of the main components of engineering structures applicable in the mechanical and electrical industries. They can be used in macro-, micro-, or nano-scales. Various analyses of the beam structures in mentioned scales need more consideration and theories. The various environments such as thermal, electrical, and different types of foundation can significantly change the responses of the beams.^1,2 In addition, manufacturing of beams with piezoelectric materials leads to intelligent systems that can be used as elements of nano-electro-mechanical-systems. Combination of the above-mentioned aspects of beams leads to an interesting problems in scope of mechanical engineering and nano-electro-mechanical-systems.^3,4 In recent years, advances in production of advanced materials have led to the emergence of new materials known as functionally graded materials (FGMs). Smooth and continuous variation of thermo-mechanical properties in a specific direction is the important feature of FGMs. For this reason, in structures made of FGMs, there are smaller stress concentrations and also reduction or elimination of cracking and delamination phenomenon.^5–12 Advances in the methods of FGM production and creation of new and intelligent structures made of FGM have encouraged various researchers for presentation of advanced analysis of these structures.^13–16

The behaviors of the beams in various subjects were studied by the researchers in the form of articles and books.^17–21 One of these main subjects is wave propagation and vibrational analysis of the beams based on various beam theories.^22,23 Various beam theories such as the Euler–Bernoulli, Timoshenko or first-order shear deformation, Reddy or parabolic shear deformation, and Levinson beam theories were employed based on the nonlocal differential constitutive relations of Eringen by Reddy.²⁴ Li et al.²⁵ developed an analytical model of small-scaled functionally graded (FG) beams for the flexural wave propagation analysis based on the nonlocal strain gradient theory (NSGT). The size-dependent wave propagation analysis of double-piezoelectric nano-beam-systems (DPNBSs) based on Euler–Bernoulli beam model was carried out by Ghorbanpour Arani et al.²⁶ They concluded that the imposed external voltage is an effective controlling parameter for wave propagation of the coupled system. Ke et al. investigated the dispersion behavior of waves in magneto-electro-elastic nano-beams based on Euler model and Timoshenko nano-beam model.²⁷ In the other research, a sandwich beam with periodic multiple dissipative resonators in the sandwich core material was investigated for broadband wave mitigation and/or absorption by Chen et al.²⁸ Joglekar and Mitra²⁹ presented an analytical–numerical method, based on the use of wavelet spectral finite elements, in order to study the nonlinear interaction of flexural waves with a breathing crack present in a slender beam. Ding et al.³⁰ studied the propagation and attenuation properties of waves in ordered and disordered periodic composite Timoshenko beams, which consider the effects of axial static load and structural damping. They assumed that beam is resting on elastic foundations and subjected to moving loads of constant amplitude with a constant velocity.

Recently, the influences of the length scales parameters were considered by the researchers in nano- and micro-structures. In the aforementioned studies in order to incorporate the small scales in equations of motions, various theories such as the strain gradient theory and Eringen’s differential nonlocal model were used.^31,32 Classical continuum models,^33,34 nonlocal continuum theory,^35–37 strain gradient theory,^38,39 and modified couple stress models^40–42 have been used by researchers for analysis of nano/micro systems. Generally, based on the nonlocal continuum theory, the stress at a specified point of the body depends on the strains at other near points.⁴³ In addition, with regard to the strain gradient theory, the strain energy is a function of the strain and curvature components tensors.⁴⁴ Lim et al.⁴⁵ carried out an investigation to show that the nonlocal and strain gradient parameters basically described two different physical properties of the structures in nano- and micro-scales. They have presented a new approach and theory to relate the strain gradient and nonlocal theories named as NSGT. This theory is a combination of the two aforementioned theories that incorporated both small scales parameters, namely, nonlocal and strain gradient parameters. Based on the NSGT, Liew et al.⁴⁶ analyzed the wave propagation in a single-walled carbon nanotubes by molecular dynamics simulations

Electromechanical coupling between electric fields and mechanical fields has been widely explored in the development of piezoelectricity-based devices in transduction technology, including transducers, sensors, resonators, and energy harvesters. In order to fulfill these applications of piezoelectric nanostructures, it is essential to get a thorough and comprehensive understanding on the electromechanical coupling behaviors of piezoelectric materials at the nano-scale. On the other hand, the flexoelectricity that refers to the spontaneous polarization in response to a nonuniform strain or a strain gradient is a universal effect in all classes of dielectric materials even in the centrosymmetric crystals. According to the aforementioned wide applications and the importance of the subject, the flexoelectricity property is one of the new and interesting subjects in recent times, and many researchers are studying on this subject.

The flexoelectricity is related to a particular electromechanical coupling phenomenon between polarization and strain gradients.⁴⁷ Zhang et al.⁴⁸ investigated the flexoelectric effect on the electro-elastic responses and the free vibrational behaviors of a piezoelectric nano-plate (PNP). Also, a modified Kirchhoff plate model considering the residual surface stress, surface elasticity, surface piezoelectricity, and flexoelectricity was developed to investigate the electro-elastic responses and vibrational behaviors of a bending PNP by Zhang and Jiang.⁴⁹ Other researchers in several investigations studied effects of the flexoelectricity on various behaviors of the beam.^50–52 On the other hand, surface effects have significant influence on nano- and micro-structures.⁵³ Nonlinear free vibration of simply supported FG nano-scale beams with considering surface effects (surface elasticity, tension, and density) and balance condition between the FG nano-beam bulk and its surfaces was investigated by Hosseini-Hashemi and Nazemnezhad.⁵⁴ Free vibration of current-carrying nano-scaled beams incorporating Rayleigh, Timoshenko, and higher order beam models considering the surface energy were investigated by Kiani.⁵⁵ He studied the effects of surface and shear deformation, electric current, magnetic field strength, and geometric parameters of the nano-beam on the first 10 natural frequencies.

With attention to the literature review mentioned above and author’s knowledge, we can conclude that there is no published work about nonlinear vibration analysis of the FG nano-beams based on the nonlocal elasticity theory considering surface and flexoelectric effects. The nano-beam is rested on Winkler–Pasternak foundation. The material properties of the nano-beam are assumed to change continuously along the thickness direction according to different patterns of material distribution. The effects of surface elasticity, di-electricity, and piezoelectricity as well as bulk flexoelectricity are included in equations to consider coupling of strain gradients and electrical polarizations. First-order shear deformation beam theory (FSDBT), residual surface stresses, and Hamilton’s principle are used to derive the governing equations of motion. The analytical method is used to calculate the nonlinear natural frequency of flexoelectric FG Timoshenko nano-beam as well as nonlinear mode shape. Furthermore, the comprehensive numerical results are presented to investigate the influence of important parameters such as flexoelectricity of the surface and bulk, residual surface stresses, Winkler and shear coefficients, power index, and geometric parameters on the nonlinear vibration behaviors of flexoelectric FG nano-beam.

Material properties of FG nano-beams

In this section, the material properties of flexoelectric FG nano-beam are expressed in detail. The schematic representation of the problem is presented in Figure 1.

Figure 1.

Flexoelectric FG nano-beam and attached coordinate system.

According to this figure, it is assumed that different patterns are used for distribution of the material properties in FG nano-beam along the thickness direction. Based on the aforementioned patterns, the volume fraction of constituent PZH-5H, $V_{p}$ is written in the following form

\begin{array}{l} V_{p} = {(\frac{2 z + h}{2 h})}^{n}, - \frac{1}{2} h \leq z \leq - \frac{1}{2} h f o r s i m p l e p o w e r l o w (P - F G M) \\ V_{p} = \frac{1}{2} {(\frac{\frac{h}{2} + z}{\frac{h}{2}})}^{n}, - \frac{1}{2} h \leq z \leq 0, V_{p} = (1 - \frac{1}{2} {(\frac{\frac{h}{2} - z}{\frac{h}{2}})}^{n}), 0 \leq z \leq \frac{1}{2} h f o r S i g m o i d l a w (S - F G M) \end{array}

(1)

In equation (1), $h$ is the thickness of flexoelectric FG nano-beam and $z$ is the coordinate along thickness direction. Also, $n$ is known as power gradient index for FG materials. Figure 2 shows the distribution of PZT-5H in which $V_{p}$ is the volume fraction along the thickness direction.

Figure 2.

Distribution of materials volume fraction along thickness coordinates for a FG nano-beam. (a) simple power low (P-FGM) and (b) Sigmoid law (S-FGM).

The nonhomogeneous material properties of flexoelectric FG nano-beam are obtained using the mixture Voigt rule.^56,57 Based on the aforementioned rule, a symbolic material property $P (z)$ is assumed as

\begin{array}{l} P (z) = P_{b} + P_{p b} {(\frac{2 z + h}{2 h})}^{n}, - \frac{1}{2} h \leq z \leq \frac{1}{2} h . s i m p l e p o w e r l o w (P - F G M) \\ P_{pb} = P_{p} - P_{b} \\ P (z) = P_{p} \exp (- \frac{1}{2} \ln (\frac{P_{p}}{P_{b}}) (1 - \frac{2 z}{h})), - \frac{1}{2} h \leq z \leq \frac{1}{2} h . E x p o n e n t i a l l a w (E - F G M) \end{array} \begin{array}{l} P (z) = \frac{1}{2} {(\frac{\frac{h}{2} + z}{\frac{h}{2}})}^{n} P_{p} + (1 - \frac{1}{2} {(\frac{\frac{h}{2} + z}{\frac{h}{2}})}^{n}) P_{b}, - \frac{1}{2} h \leq z \leq 0. \\ P (z) = (1 - \frac{1}{2} {(\frac{\frac{h}{2} - z}{\frac{h}{2}})}^{n}) P_{p} + \frac{1}{2} {(\frac{\frac{h}{2} + z}{\frac{h}{2}})}^{n} P_{b}, 0 \leq z \leq \frac{1}{2} h . S i g m o i d l a w (S - F G M) \end{array}

(2)

where

P_{p}

and

P_{b}

are the PZH-5H and

{B a T i o}_{3}

properties, respectively, and also

P_{p b} = P_{p} - P_{b}

. In the present work, it is assumed that the elasticity modulus

E

and density

ρ

are described by equation (2), while Poisson’s ratio

υ

is considered to be constant.

Formulations

For a nano-dielectric material including flexoelectricity effect, the electric Gibbs free energy density function $U_{b}$ can be written as^48,49

U_{b} = - \frac{1}{2} a_{i j} E_{i} E_{j} + \frac{1}{2} C_{i j k} ε_{i j} ε_{k l} - e_{i j k} E_{i} ε_{k l} - f_{i j k l} E_{i} η_{i k l}

(3)

where

f

is the fourth-order flexoelectricity tensor,

a

is the second-order dielectric constant tensor,

c

is the fourth-order elastic constant tensor, and

e

is the classical piezoelectric tensor. It is worth mentioning that the fifth- and sixth-order terms are neglected.^48,49 It should be noted that in the present investigation, the effects of the higher order couplings between the strain and the strain gradient and the strain gradient and the polarization gradient are neglected for simplicity. Considering equation (3), the constitutive equations for the bulk can be derived as

\begin{array}{l} σ_{i j} = \frac{\partial U_{b}}{\partial ε_{i j}} = C_{i j k l} ε_{i j} - e_{i j k} E_{i} + f_{k l i j} \frac{\partial E_{k}}{\partial x_{l}} \\ τ_{j k l} = \frac{\partial U_{b}}{\partial η_{i k l}} = - f_{i j k l} E_{i} \\ D_{i} = - \frac{\partial U_{b}}{\partial E_{i}} = a_{i j} E_{j} + e_{i j k} ε_{j k} + f_{i j k l} η_{j k l} \end{array}

(4)

where

σ_{i j}

is the Cauchy stress tensor,

D_{i}

is the electric displacement vector, and

τ_{j k l}

is the higher order stress (moment stress), respectively. On the other hand, considering the surface effects such as the surface elasticity, the residual surface stress, and the surface piezoelectricity, the surface internal energy density

U_{s}

can be specified as the following form^47,52,58

U_{s} = Γ_{α β} ε_{α β}^{s} - \frac{1}{2} α_{γ κ}^{s} E_{γ}^{s} E_{κ}^{s} + \frac{1}{2} c_{α β γ κ}^{s} ε_{α β}^{s} ε_{γ κ}^{s} - e_{κ α β}^{s} E_{κ}^{s} ε_{α β}^{s}

(5)

In above equation, $Γ_{α β}, α_{γ κ}^{s}, c_{α β γ κ}^{s}, e_{κ α β}^{s},$ and $E_{κ}$ are characterized as the surface residual stress tensor, the surface permittivity, elastic constants, surface piezoelectric tensor, and surface electric field, respectively. Considering equation (5), the constitutive equations for the surface of the beam can be derived as

\begin{array}{l} σ_{i j}^{s} = \frac{\partial U_{s}}{\partial ε_{i j}^{s}} = Γ_{α β} + C_{α β γ κ}^{s} ε_{γ κ}^{s} - e_{κ α β}^{s} E_{κ}^{s} \\ D_{γ}^{s} = - \frac{\partial U_{s}}{\partial E_{γ}^{s}} = a_{γ κ}^{s} E_{κ}^{s} + e_{γ α β}^{s} ε_{j k} \end{array}

(6)

where

σ_{α β}^{s}

and

D_{γ}^{s}

are the surface Cauchy stress and surface electric displacement. The displacement field based on the first-order shear deformation beam theory (FSDBT) is expressed as

\begin{array}{l} u (x, z, t) = u_{0} (x, t) - z ϕ (x, t) \\ w (x, z, t) = w_{0} (x, t) \end{array}

(7)

where

u

and

w

are axial and transverse displacement components and

u_{0}

and

w_{0}

are axial and transverse displacement of mid-surface. In addition,

ϕ

is the rotation angle of the beam cross-section due to the pure bending, which in Euler–Bernoulli beam is defined as

ϕ = \frac{\partial w}{\partial x}

. The relations between displacement and strain as well as strain gradient are specified as the following forms

ε_{i j} = \frac{1}{2} (u_{j, i} + u_{i, j}), η_{j k l} = ε_{k l, j} = \frac{1}{2} (u_{k, j i} + u_{l, j k})

(8)

where

u_{i}

are the components of the mechanical displacement. In addition, for surface of the beam, strain–displacement relation is defined:

ε_{α β}^{s} = \frac{1}{2} (u_{β, α}^{s} + u_{α, β}^{s})

(9)

Regarding equation (8), nonzero strains and strain gradients are calculated using equation (7) as follows

\begin{array}{l} ε_{x x} = \frac{\partial u_{0}}{\partial x} - z \frac{\partial ϕ}{\partial x}, \end{array} \begin{array}{l} ε_{x z} = \frac{1}{2} (\frac{\partial w}{\partial x} - ϕ), \\ η_{111} = \frac{\partial ε_{x x}}{\partial x} = \frac{\partial^{2} u_{0}}{\partial x^{2}} - z \frac{\partial^{2} ϕ}{\partial x^{2}}, \\ η_{113} = \frac{\partial ε_{x x}}{\partial z} = - \frac{\partial ϕ}{\partial x}, \\ η_{131} = η_{311} = \frac{\partial ε_{x z}}{\partial x} = \frac{1}{2} (\frac{\partial^{2} w}{\partial x^{2}} - \frac{\partial ϕ}{\partial x}) \end{array}

(10)

According to the NSGT,^12,59 the constitutive relations are rewritten as the following form

\begin{array}{l} (1 - {(e_{0} a)}^{2} \nabla^{2}) σ_{i j} = \frac{\partial U_{b}}{\partial ε_{i j}} = C_{i j k l} (1 - l_{m}^{2} \nabla^{2}) ε_{k l} - e_{k i j} E_{k} + \frac{1}{2} f_{k l i j} \frac{\partial E_{k}}{\partial x_{l}} \\ (1 - {(e_{0} a)}^{2} \nabla^{2}) τ_{j k l} = \frac{\partial U_{b}}{\partial η_{i k l}} = - \frac{1}{2} f_{i j k l} E_{i} \\ (1 - {(e_{0} a)}^{2} \nabla^{2}) D_{i} = - \frac{\partial U_{b}}{\partial E_{i}} = a_{i j} E_{j} + e_{i j k} (1 - l_{m}^{2} \nabla^{2}) ε_{j k} + \frac{1}{2} f_{i j k l} η_{j k l} \\ (1 - {(e_{0} a)}^{2} \nabla^{2}) σ_{i j}^{s} = \frac{\partial U_{s}}{\partial ε_{i j}^{s}} = Γ_{α β} + C_{α β γ κ}^{s} (1 - l_{m}^{2} \nabla^{2}) ε_{γ κ}^{s} - e_{κ α β}^{s} E_{κ}^{s} \\ (1 - {(e_{0} a)}^{2} \nabla^{2}) D_{γ}^{s} = - \frac{\partial U_{s}}{\partial E_{γ}^{s}} = a_{γ κ}^{s} E_{κ}^{s} + e_{γ α β}^{s} (1 - l_{m}^{2} \nabla^{2}) ε_{j k} \end{array}

(11)

In equation (11), $\nabla^{2} = \partial^{2} / \partial x^{2}$ is the Laplacian operator, $l_{m}$ is the strain gradient parameter, and $e_{0} a$ is the nonlocal parameter. The electric field is assumed along the thickness direction due to small ratio of thickness to length which means that $E_{1} = E_{2} = 0$ and $E_{3} \neq 0$ . It is worth mentioning that with the open circuit condition, the electric displacement on the surface is zero. Using equations (10) and (11), the electric field can be obtained as

E_{3} = - \frac{e_{31}}{a_{33}} \frac{\partial u_{0}}{\partial x} + (z \frac{e_{31}}{a_{33}} + \frac{1}{2} \frac{f_{111}}{a_{33}} + \frac{1}{2} \frac{f_{14}}{a_{33}}) \frac{\partial ϕ}{\partial x} - \frac{1}{2} \frac{f_{14}}{a_{33}} \frac{\partial^{2} w}{\partial x^{2}}

(12)

And consequently the electric field gradient

\frac{\partial E_{3}}{\partial z} = \frac{e_{31}}{a_{33}} \frac{\partial ϕ}{\partial x}

(13)

According to equations (10), (11), and (13), the nonzeros stress and moment stresses can be calculated as the following forms

\begin{array}{l} (1 - {(e_{0} a)}^{2} \nabla^{2}) σ_{11} = ((1 - l_{m}^{2} \nabla^{2}) C_{11} + \frac{e_{31}^{2}}{a_{33}}) \frac{\partial u}{\partial x} + (\frac{1}{2} f_{14} \frac{e_{31}}{a_{33}} - (1 - l_{m}^{2} \nabla^{2}) C_{11} z - e_{31} (z \frac{e_{31}}{a_{33}} + \frac{1}{2} \frac{f_{111}}{a_{33}} + \frac{1}{2} \frac{f_{14}}{a_{33}})) \frac{\partial ϕ}{\partial x} \\ + \frac{1}{2} e_{31} \frac{f_{14}}{a_{33}} \frac{\partial^{2} w}{\partial x^{2}} - \frac{e_{31}}{2 a_{33}} {(\frac{\partial w}{\partial x})}^{2} + (1 - l_{m}^{2} \nabla^{2}) C_{11} {(\frac{\partial w}{\partial x})}^{2} \\ (1 - {(e_{0} a)}^{2} \nabla^{2}) σ_{13} = k (1 - l_{m}^{2} \nabla^{2}) C_{44} (\frac{\partial w}{\partial x} - ϕ) + \frac{1}{2} f_{111} \frac{e_{31}}{a_{33}} \frac{\partial ϕ}{\partial x} \\ (1 - {(e_{0} a)}^{2} \nabla^{2}) τ_{113} = \frac{1}{2} \frac{e_{31} f_{14}}{a_{33}} \frac{\partial u}{\partial x} - \frac{1}{2} (\frac{f_{111} f_{14}}{a_{33}} + \frac{f_{14}^{2}}{a_{33}} + \frac{e_{31} f_{14}}{a_{33}}) \frac{\partial ϕ}{\partial x} + \frac{1}{2} \frac{f_{14}^{2}}{a_{33}} \frac{\partial^{2} w}{\partial x^{2}} + \frac{1}{2} \frac{e_{31} f_{14}}{a_{33}} {(\frac{\partial w}{\partial x})}^{2} \\ (1 - {(e_{0} a)}^{2} \nabla^{2}) τ_{131} = (1 - {(e_{0} a)}^{2} \nabla^{2}) τ_{311} = \frac{1}{2} \frac{e_{31} f_{111}}{a_{33}} \frac{\partial u}{\partial x} - \frac{1}{2} (\frac{f_{111} f_{14}}{a_{33}} + \frac{f_{111}^{2}}{a_{33}} + \frac{e_{31} f_{111}}{a_{33}}) \frac{\partial ϕ}{\partial x} \\ + \frac{1}{2} \frac{f_{111} f_{14}}{a_{33}} \frac{\partial^{2} w}{\partial x^{2}} + \frac{1}{2} \frac{e_{31} f_{111}}{a_{33}} {(\frac{\partial w}{\partial x})}^{2} \end{array}

(14)

where

k

is the shear correction factor, which in this investigation is selected

5 / 6

. The Hamilton’s principle is used to drive the governing equation of motion as follows⁶⁰

0 = \int_{0}^{T} (δ T - δ U_{s} - δ U_{f} + δ W) d \bar{t}

(15)

where

{δ U}_{s}

{δ U}_{f}, δ T

, and

δ W

are the variations of strain energy, foundation reaction, kinetic energy, and external works, respectively. Variation of strain energy

δ U_{s}

is calculated as

δ U_{s} = \int_{0}^{L} \int_{A} (σ_{11} δ ε_{11} + 2 σ_{13} δ ε_{13} + τ_{113} δ η_{113} + 2 τ_{131} δ η_{131}) d A d x

(16)

Variation of kinetic energy is represented as

δ T = \int_{0}^{L} \int_{A} ρ (z) (\frac{\partial u}{\partial t} \frac{\partial u}{\partial t} + \frac{\partial w}{\partial t} \frac{\partial w}{\partial t}) d A d x

(17)

Variations of work done by the external forces and Winkler–Pasternak foundation are written as^5,60,61

\begin{array}{l} δ W = \int_{0}^{L} (F δ u_{0} + Q δ w_{0} + {\bar{N}}_{0} \frac{\partial w_{0}}{\partial \bar{x}} \frac{\partial δ w_{0}}{\partial \bar{x}}) d x \\ δ U_{f} = \int_{0}^{L} \int_{0}^{b} (K_{w} w_{0} δ w_{0} + K_{s} \frac{\partial w_{0}}{\partial x} δ (\frac{\partial w_{0}}{\partial x}) + K_{n l} w_{0}^{3} δ w_{0}) d \bar{y} d \bar{x} \end{array}

(18)

where F and Q are the axial and transverse forces per unit length, respectively, and

{\bar{N}}_{0}

is the axial compressive or pretension force. Also,

K_{w}

K_{n l,}

and

K_{s}

are linear stiffness, nonlinear stiffness, and shear coefficients of the foundation, respectively. Substituting equations (10) and (14) into equations (16) and (18) and consequently into equation (15) yields the governing equations of motions as follows

\begin{array}{l} δ u_{0} : \frac{\partial N_{x x}}{\partial x} + \frac{\partial N_{x x}^{s}}{\partial x} + F = I_{0} \frac{\partial^{2} u_{0}}{\partial t^{2}} - I_{1} \frac{\partial^{2} ϕ}{\partial t^{2}} \\ δ w_{0} : \frac{\partial N_{x z}}{\partial x} - \frac{\partial^{2} {\bar{N}}_{131}}{\partial x^{2}} - {\bar{K}}_{w} w_{0} + {\bar{K}}_{s} \frac{\partial^{2} w_{0}}{\partial x^{2}} - {\bar{K}}_{n l} w_{0}^{3} + σ_{0}^{s} \frac{\partial^{2} w_{0}}{\partial x^{2}} + Q = I_{0} \frac{\partial^{2} w_{0}}{\partial t^{2}} \\ δ ϕ : - \frac{\partial M_{x x}^{s}}{\partial x} + N_{x z} - \frac{\partial M_{x x}}{\partial x} - \frac{\partial {\bar{N}}_{131}}{\partial x} - \frac{\partial {\bar{N}}_{113}}{\partial x} = I_{1} \frac{\partial^{2} u_{0}}{\partial t^{2}} + I_{2} \frac{\partial^{2} ϕ}{\partial t^{2}} \end{array}

(19)

where

N_{\bar{x}}

Q_{\bar{x z}}

M_{\bar{x}}

, and

M_{\bar{x}}^{h}

are the resultants of forces and moments and

σ_{0}

is the residual surface stress.⁴⁷ The nonlocal resultants of forces and the moments are expressed as

\begin{array}{l} (1 - {(e_{0} a)}^{2} \nabla^{2}) N_{x x} = A_{x 1} \frac{\partial u}{\partial x} - A_{x 2} \frac{\partial^{3} u}{\partial x^{3}} + A_{x 3} \frac{\partial ϕ}{\partial x} - B_{x 1} \frac{\partial ϕ}{\partial x} + B_{x 2} \frac{\partial^{3} ϕ}{\partial x^{3}} + A_{x 4} \frac{\partial^{2} w}{\partial x^{2}} + A_{x 5} {(\frac{\partial w}{\partial x})}^{2} \\ - \frac{A_{x 2}}{2} ({(\frac{\partial^{2} w}{\partial x^{2}})}^{2} + \frac{\partial w}{\partial x} \frac{\partial^{3} w}{\partial x^{3}}), \\ (1 - {(e_{0} a)}^{2} \nabla^{2}) M_{x x} = B_{x 1} \frac{\partial u}{\partial x} - B_{x 2} \frac{\partial^{3} u}{\partial x^{3}} + B_{x 3} \frac{\partial ϕ}{\partial x} - D_{x 1} \frac{\partial ϕ}{\partial x} + D_{x 2} \frac{\partial^{3} ϕ}{\partial x^{3}} + B_{x 4} \frac{\partial^{2} w}{\partial x^{2}} + B_{x 5} {(\frac{\partial w}{\partial x})}^{2} \\ - \frac{B_{x 2}}{2} ({(\frac{\partial^{2} w}{\partial x^{2}})}^{2} + \frac{\partial w}{\partial x} \frac{\partial^{3} w}{\partial x^{3}}), \\ (1 - {(e_{0} a)}^{2} \nabla^{2}) N_{x z} = A_{x z 1} (\frac{\partial w}{\partial x} - ϕ) + A_{x z 2} \frac{\partial ϕ}{\partial x} - A_{x z 3} (\frac{\partial^{3} w}{\partial x^{3}} - \frac{\partial ϕ}{\partial x}), \end{array} \begin{array}{l} (1 - {(e_{0} a)}^{2} \nabla^{2}) {\bar{N}}_{131} = A_{131}^{1} \frac{\partial u}{\partial x} - (A_{131}^{2} + B_{131}^{1}) \frac{\partial ϕ}{\partial x} + A_{131}^{3} \frac{\partial^{2} w}{\partial x^{2}} + A_{131}^{4} {(\frac{\partial w}{\partial x})}^{2}, \\ (1 - {(e_{0} a)}^{2} \nabla^{2}) {\bar{N}}_{113} = A_{113}^{1} \frac{\partial u}{\partial x} - (A_{113}^{2} + B_{113}^{1}) \frac{\partial ϕ}{\partial x} + A_{113}^{3} \frac{\partial^{2} w}{\partial x^{2}} + A_{113}^{4} {(\frac{\partial w}{\partial x})}^{2}, \\ (1 - {(e_{0} a)}^{2} \nabla^{2}) N_{x x}^{s} = A^{s}_{x 1} \frac{\partial u}{\partial x} - A^{s} x_{2} \frac{\partial^{3} u}{\partial x^{3}} + A^{s}_{x 3} \frac{\partial ϕ}{\partial x} - B^{s}_{x 1} \frac{\partial ϕ}{\partial x} + B^{s} x_{2} \frac{\partial^{3} ϕ}{\partial x^{3}} + A^{s}_{x 4} \frac{\partial^{2} w}{\partial x^{2}} + A^{s} x_{5} {(\frac{\partial w}{\partial x})}^{2} \\ - \frac{A^{s} x_{2}}{2} ({(\frac{\partial^{2} w}{\partial x^{2}})}^{2} + \frac{\partial w}{\partial x} \frac{\partial^{3} w}{\partial x^{3}}), \\ (1 - {(e_{0} a)}^{2} \nabla^{2}) M_{x x}^{s} = B^{s}_{x 1} \frac{\partial u}{\partial x} - B^{s} x_{2} \frac{\partial^{3} u}{\partial x^{3}} + B^{s}_{x 3} \frac{\partial ϕ}{\partial x} - D^{s}_{x 1} \frac{\partial ϕ}{\partial x} + D^{s} x_{2} \frac{\partial^{3} ϕ}{\partial x^{3}} + B^{s}_{x 4} \frac{\partial^{2} w}{\partial x^{2}} + B^{s} x_{5} {(\frac{\partial w}{\partial x})}^{2} \\ - \frac{B^{s} x_{2}}{2} ({(\frac{\partial^{2} w}{\partial x^{2}})}^{2} + \frac{\partial w}{\partial x} \frac{\partial^{3} w}{\partial x^{3}}) \end{array}

(20)

It is noted that, $A_{x}$ , $B_{x,}$ and $D_{x}$ in equation (20) are the stretching stiffness, stretching-bending coupling stiffness, and bending stiffness coefficients, respectively, which can be found in Appendix. The integration of other constants presented in equation (20) can be calculated as

{\begin{matrix} I_{0} \\ I_{1} \\ I_{2} \end{matrix}} = \int_{- h / 2}^{h / 2} {\begin{matrix} 1 \\ z \\ z^{2} \end{matrix}} ρ (z) d z

(21)

To obtain the equations of motion, equation (20) should be substituted into equation (19). Therefore, three coupled equations of motion are obtained

\begin{array}{l} δ u_{0} : (A_{x 1} + A_{x 1}^{s}) \frac{\partial^{2} u_{0}}{\partial x^{2}} - (A_{x 2} + A_{x 2}^{s}) \frac{\partial^{4} u_{0}}{\partial x^{4}} - (B_{x 1} + B_{x 1}^{s}) \frac{\partial^{2} ϕ}{\partial x^{2}} + (B_{x 2} + B_{x 2}^{s}) \frac{\partial^{4} ϕ}{\partial x^{4}} + (A_{x 3} + A_{x 3}^{s}) \frac{\partial^{2} ϕ}{\partial x^{2}} \\ + (A_{x 4} + A_{x 4}^{s}) \frac{\partial^{3} w}{\partial x^{3}} + (A_{x 5} + A_{x 5}^{s}) {(\frac{\partial w}{\partial x})}^{2} - (A_{x 2} + A_{x 2}^{s}) ({(\frac{\partial^{2} w}{\partial x^{2}})}^{2} + \frac{\partial w}{\partial x} \frac{\partial^{3} w}{\partial x^{3}}) \\ = (1 - {(e_{0} a)}^{2} \nabla^{2}) I_{0} \frac{\partial^{2} u_{0}}{\partial t^{2}} - (1 - {(e_{0} a)}^{2} \nabla^{2}) I_{1} \frac{\partial^{2} ϕ}{\partial t^{2}} \\ δ w_{0} : - A_{131}^{1} \frac{\partial^{3} u_{0}}{\partial x^{3}} - A_{x z 1} \frac{\partial ϕ}{\partial x} - A_{x z 3} (\frac{\partial^{4} w}{\partial x^{4}} - \frac{\partial^{2} ϕ}{\partial x^{2}}) + A_{x z 2} \frac{\partial^{2} ϕ}{\partial x^{2}} + (B_{131}^{1} + A_{131}^{2}) \frac{\partial^{3} ϕ}{\partial x^{3}} + A_{x z 1} \frac{\partial^{2} w}{\partial x^{2}} - A_{131}^{3} \frac{\partial^{4} w}{\partial x^{4}} \\ - A_{131}^{4} ({(\frac{\partial^{2} w}{\partial x^{2}})}^{2} + \frac{\partial w}{\partial x} \frac{\partial^{3} w}{\partial x^{3}}) + \frac{\partial^{2} w}{\partial x^{2}} ((A_{x 1} + A_{x 1}^{s}) \frac{\partial^{2} u_{0}}{\partial x^{2}} - (A_{x 2} + A_{x 2}^{s}) \frac{\partial^{4} u_{0}}{\partial x^{4}} - (B_{x 1} + B_{x 1}^{s}) \frac{\partial^{2} ϕ}{\partial x^{2}} \\ + (B_{x 2} + B_{x 2}^{s}) \frac{\partial^{4} ϕ}{\partial x^{4}} + (A_{x 3} + A_{x 3}^{s}) \frac{\partial^{2} ϕ}{\partial x^{2}} + (A_{x 4} + A_{x 4}^{s}) \frac{\partial^{3} w}{\partial x^{3}} + (A_{x 5} + A_{x 5}^{s}) {(\frac{\partial w}{\partial x})}^{2} \\ - (A_{x 2} + A_{x 2}^{s}) ({(\frac{\partial^{2} w}{\partial x^{2}})}^{2} + \frac{\partial w}{\partial x} \frac{\partial^{3} w}{\partial x^{3}})) = (1 - {(e_{0} a)}^{2} \nabla^{2}) I_{0} \frac{\partial^{2} w}{\partial t^{2}} \\ δ ϕ : (- B_{x 2} - B_{x 2}^{s} - A_{131}^{1} - A_{113}^{1}) \frac{\partial^{2} u_{0}}{\partial x^{2}} - (B_{x 2} + B_{x 2}^{s}) \frac{\partial^{4} u_{0}}{\partial x^{4}} - A_{x z 1} ϕ + A_{x z 2} \frac{\partial ϕ}{\partial x} + + (D_{x 2} + D_{x 2}^{s}) \frac{\partial^{4} ϕ}{\partial x^{4}} \\ (- B_{x 3} + D_{x 1} + D_{x 1}^{s} - B_{x 3}^{s} + B_{131}^{1} + A_{131}^{1} + B_{113}^{1} + A_{113}^{1}) \frac{\partial^{2} ϕ}{\partial x^{2}} + A_{x z 1} \frac{\partial w}{\partial x} - A_{x z 1} ϕ \\ - (A_{131}^{4} + A_{113}^{4} + B_{x 5} + B_{x 5}^{s}) \frac{d}{d x} ({(\frac{\partial w}{\partial x})}^{2}) - (B_{x 4} + B_{x 4}^{s} + A_{131}^{3} + A_{113}^{3}) \frac{\partial^{3} w}{\partial x^{3}} + A_{x z 2} \frac{\partial ϕ}{\partial x} \\ + (B_{x 2} + B_{x 2}^{s}) ({(\frac{\partial^{2} w}{\partial x^{2}})}^{2} + \frac{\partial w}{\partial x} \frac{\partial^{3} w}{\partial x^{3}})) = - (1 - {(e_{0} a)}^{2} \nabla^{2}) I_{0} \frac{\partial^{2} u_{0}}{\partial t^{2}} + (1 - {(e_{0} a)}^{2} \nabla^{2}) I_{2} \frac{\partial^{2} ϕ}{\partial t^{2}} \end{array}

(22)

Boundary conditions for solution of equation (31) are defined as⁶²

\begin{array}{l} u = w = ϕ = Φ = 0 f o r c l a m p e d e n d s \\ u = w = \frac{\partial ϕ}{\partial x} = Φ = 0 f o r h i n g e d e n d s \end{array}

(23)

Nonlinear vibration analysis

In order to solve the nonlinear governing equations of motion, differential quadrature (DQ) method is employed. Based on the aforementioned method, the approximate solution of a function $f (x)$ can be found in the following form

f (x) = \sum_{j = 1}^{N} λ_{j} ψ_{j} (x)

(24)

where

N

is the total number of grid points inside a closed interval. In this method, the smooth basis functions are selected as the various functions form such as Chebyshev polynomials, exponential polynomials, and Fourier polynomials.^63,64 Also, equation (24) for the one-dimensional case can be written in matrix form as

f = C λ

(25)

where

f = {[f (x_{1}), f (x_{2}), f (x_{3}), \dots, f (x_{N})]}^{T}

is the vector of the unknown function values,

λ

is the vector of the unknown coefficients, and

λ_{j}

are the components of the coefficient matrix

C

are given by

C_{i j} = ψ_{j} (x_{i})

for

i, j = 1.2 . 3, \dots, N

.⁶⁴ It should be noted that

λ

is the independent variable that is graded in its direction by numbering

N_{λ}

. In fact,

N_{λ}

is the number of nodes in a given direction such as

λ

. The n'th order derivative of equation (24) can be computed as⁶³

\frac{d^{n} f (x)}{d x^{n}} = \sum_{j = 1}^{N} λ_{j} \frac{d^{n} ψ_{j} (x)}{d x^{n}}, f o r n = 1, 2, 3, …, N - 1

(26)

Equation (25) is rewritten as the following matrix form

f^{(n)} = C^{(n)} λ w i t h C_{i j}^{(n)} = {\frac{d^{n} ψ_{j} (x)}{d x^{n}} |}_{x_{i}} = ψ_{j}^{(n)} (x_{i}) f o r i, j = 1, 2, 3, …, N

(27)

Therefore, the governing equations and boundary conditions are discretized by means of the aforementioned method.^61,65 In this investigation, the cosine pattern is employed to generate the DQ point system as the following form

x_{j} = \frac{1}{2} {1 - \cos (\frac{π (j - 1)}{N - 1})}, j = 1, 2, …, N

(28)

In addition, column vectors for variables $u$ , $w$ , $ϕ$ , and $Φ$ are considered as follows

\begin{array}{l} u = [\begin{matrix} u_{1} & u_{2} & … & u_{N} \end{matrix}], w = [\begin{matrix} w_{1} & w_{2} & … & w_{N} \end{matrix}], \\ ϕ = [\begin{matrix} ϕ_{1} & ϕ_{2} & … & ϕ_{N} \end{matrix}] \end{array}

(29)

To solve equation (22) and associated boundary conditions (equation (23)) for nonlinear vibration analysis of the sandwich nano-beam by DQ method, the weighting coefficients for the second, third, and fourth derivatives with attention to equation (26) are determined as the following form

{\frac{\partial^{r} f (ς, λ)}{\partial ς^{r}} |}_{(ς, λ) = (ς_{i}, λ_{j})} = \sum_{j = 1}^{N_{ς}} C_{i j}^{ς (r)} f (ς_{j}, λ_{m}) = \sum_{j = 1}^{N_{ς}} C_{i j}^{ς (r)} f_{j m}, {\begin{cases} i = 1, 2, …, N_{ς} \\ m = 1, 2, …, N_{λ} \\ r = 1, 2, …, N_{ς} - 1 \end{cases}

(30)

In which weighting coefficients $C_{i j}^{ζ}$ are expressed as

\begin{matrix} {\frac{\partial^{r} f (ς, λ)}{\partial ς^{r}} |}_{(ς, λ) = (ς_{i}, λ_{j})} = \sum_{j = 1}^{N_{ς}} C_{i j}^{ς (r)} f (ς_{j}, λ_{m}) = \sum_{j = 1}^{N_{ς}} C_{i j}^{ς (r)} f_{j m}, {\begin{cases} i = 1, 2, …, N_{ς} \\ m = 1, 2, …, N_{λ} \\ r = 1, 2, …, N_{ς} - 1 \end{cases} \\ C_{i j}^{ς} = {\begin{cases} \frac{M (ς_{i})}{(ς_{i} - ς_{j}) M (ς_{j})} f o r i \neq j \\ - \sum_{_{i \neq j}^{j = 1}}^{N_{ς}} C_{i j}^{ς} f o r i = j \end{cases} \end{matrix}

(31)

It is noted that in equation (31), $M (ς_{i})$ is represented as the following form

M (ς_{i}) = \prod_{_{i \neq j}^{j = 1}}^{N_{ς}} (ς_{i} - ς_{j})

(32)

The weighting coefficients for various derivatives such as the second, third, and fourth derivatives are defined as

\begin{array}{l} M (ς_{i}) = \prod_{_{i \neq j}^{j = 1}}^{N_{ς}} (ς_{i} - ς_{j}) \\ C_{i j}^{ς (2)} = \sum_{k = 1}^{N ς} C_{i k}^{ς (1)} C_{k j}^{ς (1)}, \\ C_{i j}^{ς (3)} = \sum_{k = 1}^{N ς} C_{i k}^{ς (1)} C_{k j}^{ς (2)} = \sum_{k = 1}^{N ς} C_{i k}^{ς (2)} C_{k j}^{ς (1)}, \\ C_{i j}^{ς (4)} = \sum_{k = 1}^{N ς} C_{i k}^{ς (1)} C_{k j}^{ς (3)} = \sum_{k = 1}^{N ς} C_{i k}^{ς (3)} C_{k j}^{ς (1)} \end{array}

(33)

By applying equations (31) and (32) to equation (22), one can obtain a set of ordinary nonlinear equations as

\begin{array}{l} δ u : δ u_{0} : (A_{x 1} + A_{x 1}^{s}) \sum_{m = 1}^{N} C_{i m}^{(2)} u_{m} - (A_{x 2} + A_{x 2}^{s}) \sum_{m = 1}^{N} C_{i m}^{(4)} u_{m} - (B_{x 1} + B_{x 1}^{s} + A_{x 3} + A_{x 3}^{s}) \sum_{m = 1}^{N} C_{i m}^{(2)} ϕ_{m} \\ + (B_{x 2} + B_{x 2}^{s}) \sum_{m = 1}^{N} C_{i m}^{(4)} ϕ_{m} + (A_{x 4} + A_{x 4}^{s}) \sum_{m = 1}^{N} C_{i m}^{(3)} w_{m} + (A_{x 5} + A_{x 5}^{s}) \sum_{m = 1}^{N} C_{i m}^{(1)} w_{m} \sum_{m = 1}^{N} C_{i m}^{(1)} w_{m} \\ - (A_{x 2} + A_{x 2}^{s}) (\sum_{m = 1}^{N} C_{i m}^{(2)} w_{m} \sum_{m = 1}^{N} C_{i m}^{(2)} w_{m} + \sum_{m = 1}^{N} C_{i m}^{(1)} w_{m} \sum_{m = 1}^{N} C_{i m}^{(3)} w_{m}) = I_{0} \frac{\partial^{2} u}{\partial t^{2}} - I_{0} {(e_{0} a)}^{2} \sum_{m = 1}^{N} C_{i m}^{(2)} \frac{\partial^{2} u}{\partial t^{2}} \\ + I_{1} \frac{\partial^{2} ϕ}{\partial t^{2}} - I_{1} {(e_{0} a)}^{2} \sum_{m = 1}^{N} C_{i m}^{(2)} \frac{\partial^{2} ϕ}{\partial t^{2}} \end{array} \begin{array}{l} δ w_{0} : - A_{131}^{1} \sum_{m = 1}^{N} C_{i m}^{(3)} u_{m} - A_{x z 1} \sum_{m = 1}^{N} C_{i m}^{(1)} ϕ_{m} - A_{x z 3} (\sum_{m = 1}^{N} C_{i m}^{(4)} w_{m} - \sum_{m = 1}^{N} C_{i m}^{(2)} ϕ_{m}) + A_{x z 2} \sum_{m = 1}^{N} C_{i m}^{(2)} ϕ_{m} \\ + (B_{131}^{1} + A_{131}^{2}) \sum_{m = 1}^{N} C_{i m}^{(3)} ϕ_{m} + A_{x z 1} \sum_{m = 1}^{N} C_{i m}^{(2)} w_{m} - A_{131}^{3} \sum_{m = 1}^{N} C_{i m}^{(4)} w_{m} - A_{131}^{4} \sum_{m = 1}^{N} C_{i m}^{(2)} w_{m} \sum_{m = 1}^{N} C_{i m}^{(2)} w_{m} \\ - A_{131}^{4} \sum_{m = 1}^{N} C_{i m}^{(1)} w_{m} \sum_{m = 1}^{N} C_{i m}^{(3)} w_{m} + \sum_{m = 1}^{N} C_{i m}^{(2)} w_{m} ((A_{x 1} + A_{x 1}^{s}) \sum_{m = 1}^{N} C_{i m}^{(2)} u_{m} - (A_{x 2} + A_{x 2}^{s}) \sum_{m = 1}^{N} C_{i m}^{(4)} u_{m} \\ - (B_{x 1} + B_{x 1}^{s}) \sum_{m = 1}^{N} C_{i m}^{(2)} ϕ_{m} + (B_{x 2} + B_{x 2}^{s}) \sum_{m = 1}^{N} C_{i m}^{(4)} ϕ_{m} + (A_{x 3} + A_{x 3}^{s}) \sum_{m = 1}^{N} C_{i m}^{(2)} ϕ_{m} + (A_{x 4} + A_{x 4}^{s}) \sum_{m = 1}^{N} C_{i m}^{(3)} w_{m} \\ + (A_{x 5} + A_{x 5}^{s}) \sum_{m = 1}^{N} C_{i m}^{(1)} w_{m} \sum_{m = 1}^{N} C_{i m}^{(1)} w_{m} - (A_{x 2} + A_{x 2}^{s}) (\sum_{m = 1}^{N} C_{i m}^{(2)} w_{m} \sum_{m = 1}^{N} C_{i m}^{(2)} w_{m} + \sum_{m = 1}^{N} C_{i m}^{(1)} w_{m} \sum_{m = 1}^{N} C_{i m}^{(3)} w_{m})) \\ - K_{w} w_{m} + K_{s} \sum_{m = 1}^{N} C_{i m}^{(2)} w_{m} - K_{n l} w_{m}^{3} = I_{0} \frac{\partial^{2} w}{\partial t^{2}} - I_{0} {(e_{0} a)}^{2} \sum_{m = 1}^{N} C_{i m}^{(2)} \frac{\partial^{2} w_{m}}{\partial t^{2}} \\ δ ϕ : (- B_{x 2} - B_{x 2}^{s} - A_{131}^{1} - A_{113}^{1}) \sum_{m = 1}^{N} C_{i m}^{(2)} u_{m} - (B_{x 2} + B_{x 2}^{s}) \sum_{m = 1}^{N} C_{i m}^{(4)} u_{m} - 2 A_{x z 1} ϕ + A_{x z 2} \sum_{m = 1}^{N} C_{i m}^{(1)} ϕ_{m} \\ + (D_{x 2} + D_{x 2}^{s}) \sum_{m = 1}^{N} C_{i m}^{(4)} ϕ_{m} (- B_{x 3} + D_{x 1} + D_{x 1}^{s} - B_{x 3}^{s} + B_{131}^{1} + A_{131}^{1} + B_{113}^{1} + A_{113}^{1}) \sum_{m = 1}^{N} C_{i m}^{(2)} ϕ_{m} \\ + A_{x z 1} \sum_{m = 1}^{N} C_{i m}^{(1)} w_{m} - (A_{131}^{4} + A_{113}^{4} + B_{x 5} + B_{x 5}^{s}) \sum_{m = 1}^{N} C_{i m}^{(1)} w_{m} \sum_{m = 1}^{N} C_{i m}^{(2)} w_{m} - (B_{x 4} + B_{x 4}^{s} + A_{131}^{3} + A_{113}^{3}) \sum_{m = 1}^{N} C_{i m}^{(3)} w_{m} \\ + A_{x z 2} \sum_{m = 1}^{N} C_{i m}^{(1)} ϕ_{m} + (B_{x 2} + B_{x 2}^{s}) (\sum_{m = 1}^{N} C_{i m}^{(2)} w_{m} \sum_{m = 1}^{N} C_{i m}^{(2)} w_{m} + \sum_{m = 1}^{N} C_{i m}^{(1)} w_{m} \sum_{m = 1}^{N} C_{i m}^{(3)} w_{m})) = \\ - I_{0} \frac{\partial^{2} u}{\partial t^{2}} + I_{0} {(e_{0} a)}^{2} \sum_{m = 1}^{N} C_{i m}^{(2)} \frac{\partial^{2} u}{\partial t^{2}} + I_{2} \frac{\partial^{2} ϕ}{\partial t^{2}} - I_{2} {(e_{0} a)}^{2} \sum_{m = 1}^{N} C_{i m}^{(2)} \frac{\partial^{2} ϕ}{\partial t^{2}} \end{array}

(34)

The boundary conditions of the flexoelectric FG nano-beam using DQ method are expressed as

\begin{array}{l} {\begin{cases} u_{1} = w_{1} = ϕ_{1} = Φ_{1} = 0 \\ u_{N} = w_{N} = ϕ_{N} = Φ_{N} = 0 \end{cases} f o r c l a m p e d e n d s \\ {\begin{cases} u_{1} = w_{1} = \sum_{m = 1}^{N} C_{1 m}^{(1)} ϕ_{m} = Φ_{1} = 0 \\ u_{N} = w_{N} = \sum_{m = 1}^{N} C_{N m}^{(1)} ϕ_{m} = Φ_{N} = 0 \end{cases} f o r h i n g e d e n d s \end{array}

(35)

The discretized forms of the governing equations can be expressed as

M \ddot{X} + {\tilde{K}}_{l} X + {\tilde{K}}_{n l} (X) = 0

(36)

In this stage, we define the unknown dynamic displacement vector as the following form

X = {{[u_{i}]}^{T}, {[w_{i}]}^{T}, {[ϕ_{i}]}^{T}}, i = 1, 2, …, N

(37)

To investigate the dynamic problem and free vibration responses, the following time-dependent fields are assumed as

X = \bar{X} \exp (i ω t)

(38)

where

\bar{X} = \{[{\bar{u}}_{i}], [{\bar{w}}_{i}], [{\bar{ϕ}}_{i}]\}

are the vibration mode shape vector and

ω

is the nondimensional, nonlinear natural frequency. Substitution of equation (38) into equation (36) yields the nonlinear eigenvalue equation

({\tilde{K}}_{l} + {\tilde{K}}_{n l}) \bar{X} - ω^{2} M \bar{X} = 0

(39)

In this study, in order to solve equation (39), a direct iterative process is employed as follows

By ignoring the nonlinear stiffness ${\tilde{K}}_{n l}$ , the linear eigenvalue problem is solved. Solution of the linear problem leads to linear natural frequency ( $ω_{l}$ ) and the associated linear mode shape vector. The mode shape is then appropriately scaled up such that the maximum transverse displacement is equal to a given vibration amplitude $w_{m a x}$ .

Using the linear mode shape obtained in previous stage ${, \tilde{K}}_{n l}$ is calculated and new eigenvalue and the corresponding eigenvector are calculated from the updated eigenvalue system (equation (39)). This would give the nonlinear eigenvalue ( $ω_{n l}$ ) and the new eigenvector.

The eigenvector is scaled up again, and step b is repeated until the following relation is satisfied by the frequency values from the two subsequent iterations

\frac{| ω^{r + 1} - ω^{r} |}{ω^{r}} < ε_{0}

(40)

where $ε_{0}$ is a small value number and in the present investigation it is taken to be 0.1%.

Numerical results and discussion

In this section, a parametric study is implemented to indicate the influences of nonlocal parameter, power gradient index, geometric dimensions of the beam, surface effects, flexoelectric property and foundation property, and other important parameters on designing and controlling the nonlinear vibrational behaviors. The material properties and geometrical specifications of the flexoelectric FG nano-beam are presented in Table 1.

Table 1.

The material and geometrical properties of the constituent material of the flexoelectric FG nano-beam.^66,67

Material	$C_{11}$ (Pa)	$C_{44}$ (Pa)	Density (kg/m³)	$e_{31}$ C/m²	Geometric dimensions (nm)		$a_{33} (\frac{C}{V m})$
PZH-5Z	102E (9)	35.5 E (9)	7600	4.4	h	b	1.76E (−8)
${B a T i o}_{3}$	167.55E (9)	44.7 E (9)	6020	17.05	2	2h	0.79E (−8)

Validation of results

To justify the accuracy of the governing equations and corresponding numerical results, a comparison with existing reference using semi-analytical methods is presented based on Euler–Bernoulli nano-beam theory. Figure 3 shows comparison between the obtained results by solving the governing equations extracted in this study and equation of motion obtained in Ebrahimi and Barati.⁴⁷ It is worth noting that the flexoelectricity and surface effects were considered in the study by Ebrahimi and Barati.⁴⁷ According to this comparison, it is deduced that the present results are in a good agreement with the obtained results by Ebrahimi and Barati.⁴⁷

Figure 3.

Fundamental frequency versus dimensionless nonlocal parameter (e₀a).

Figures 4 and 5 illustrate fundamental linear and nonlinear mode shape of the flexoelectric FG nano-beam for various boundary conditions. According to these figures, it can be concluded that employing the nonlinear analysis leads to modification of results and obtaining more accurate mode shapes rather than linear one. In addition, it is observed that employing the nonlinear analysis gives a little more displacement rather than linear one.

Figure 4.

Fundamental nonlinear mode shape of the flexoelectric FG nano-beam with C–H boundary condition.

Figure 5.

Fundamental nonlinear mode shape of the flexoelectric FG nano-beam with H–H boundary condition.

The nonlinear natural frequency of the flexoelectric FG nano-beam for various boundary conditions, nonlocal parameters ( $e_{0} a$ ,) and initial condition ( $w_{m a x}$ ) with and without surface effects are presented in Tables 2 and 3. According to these results, it is concluded that increasing the nonlocal parameter leads to the decrease of nonlinear natural frequency due to the decrease of the bending stiffness of flexoelectric nano-beam. In addition, it can be seen that the nonlinear natural frequencies are enhanced by increasing the initial condition. Also it is worth noting that considering the surface effects cause to increase the nonlinear natural frequency. Therefore, in order to have a precise nonlinear vibrational analysis of the flexoelectric FG nano-beam, considering the surface effects is very significant.

Table 2.

The nonlinear natural frequency for various initial conditions with and without surface effects.

C–C				H–C
$W_{m a x}$	With surface effects	Without surface effects	Diff (%)	$W_{m a x}$	With surface effects	Without surface effects	Diff (%)
0.01	0.5388	0.5032	6.61	0.01	0.385	0.3618	6.03
0.02	0.5421	0.5068	6.51	0.02	0.3904	0.3675	5.87
0.03	0.5476	0.5126	6.39	0.03	0.3991	0.3766	5.64
0.04	0.5552	0.5205	6.25	0.04	0.4106	0.3885	5.38
0.05	0.5647	0.5304	6.07	0.05	0.4244	0.4027	5.11
0.06	0.5759	0.542	5.89	0.06	0.44	0.4188	4.82
0.07	0.5887	0.5552	5.69	0.07	0.4572	0.4363	4.57
0.08	0.6028	0.5697	5.49	0.08	0.4756	0.4551	4.31
0.09	0.6182	0.5854	5.31	0.09	0.495	0.4747	4.10
0.1	0.6346	0.6021	5.12	0.1	0.5152	0.4951	3.90

Table 3.

The nonlinear natural frequency for various nonlocal parameters with and without surface effects.

C–C
$e_{0} a$	With surface effects	Without surface effects	diff (%)
0.03	0.5856	0.5512	5.87
0.06	0.5766	0.5428	5.86
0.09	0.5624	0.5296	5.83
0.12	0.544	0.5125	5.79
0.15	0.5228	0.4927	5.76
0.18	0.4999	0.4713	5.72
0.21	0.4763	0.4491	5.71
0.24	0.4527	0.427	5.68
0.27	0.4296	0.4054	5.63
0.3	0.4076	0.3847	5.62

With regard to the obtained results presented in Table 4, the nonlinear natural frequency of the flexoelectric FG nano-beam is significantly increased by disregarding the flexoelectricity. Therefore, it can be concluded that considering the flexoelectricity leads to the decrease of the bending stiffness of the flexoelectric FG nano-beam. Differences between nonlinear natural frequencies obtained with and without flexoelectricity indicate that this phenomenon has significant effects on nonlinear vibrational behavior of the FG nano-beams.

Table 4.

The nonlinear natural frequency for various nonlocal parameters and boundary conditions with and without flexoelectricity.

C–C				H–C
$e_{0} a$	Without flexoelectricity	With flexoelectricity	diff (%)	$e_{0} a$	Without flexoelectricity	With flexoelectricity	diff (%)
0.01	0.3869	0.3404	12.02	0.01	0.2876	0.2573	10.54
0.02	0.3862	0.3397	12.04	0.02	0.2871	0.2568	10.55
0.03	0.385	0.3383	12.13	0.03	0.2863	0.2558	10.65
0.04	0.3835	0.3362	12.33	0.04	0.2852	0.2547	10.69
0.05	0.3815	0.3337	12.53	0.05	0.2838	0.2533	10.75
0.06	0.379	0.3307	12.74	0.06	0.2821	0.2514	10.88
0.07	0.3762	0.3274	12.97	0.07	0.2802	0.2495	10.96
0.08	0.3731	0.3236	13.27	0.08	0.278	0.2471	11.12
0.09	0.3696	0.3195	13.56	0.09	0.2755	0.2445	11.25
0.1	0.3658	0.3149	13.91	0.1	0.2729	0.2418	11.40

In order to investigate the effects of different patterns for materials distribution on nonlinear vibrational behavior of the flexoelectric FG nano-beams, nonlinear natural frequencies are calculated in Table 5 for three distribution patterns of materials. With regard to the results obtained by different distribution patterns, it is noted that exponential FGM calculated the smaller values for nonlinear natural frequency of the flexoelectric FG nano-beam with respect to simple power and sigmoid patterns in various boundary conditions. In addition, the nonlinear natural frequencies computed by simple power and sigmoid FGMs distribution approximately have near values. Among various boundary conditions considered in this study, C–C has the larger values with respect to other boundary condition such as H–C and S–S.

Table 5.

The nonlinear natural frequency for various FGM distributions.