Electro-thermal buckling of FG graphene platelets-strengthened piezoelectric beams under humid conditions

Abstract

In this article, thermal buckling analysis of multilayer functionally graded graphene platelets (FGGPL) strengthened piezoelectric beam subjected to external electric voltage as well as humid conditions is illustrated. The effective Young’s modulus of the nanocomposite beam is estimated within the framework of Halpin-Tsai model. While, Poisson’s ratio, mass density, and piezoelectric properties are calculated by the rule of the mixture. Four FGGPL distribution types are considered in this study. A refined two-unknown beam theory considering shear deformation as well as thickness stretching effect is employed to describe the displacement components. The principle of virtual work including thermal, moisture, and electric loads is used to derive the stability differential equations. To check the accuracy of the obtained buckling temperature, some comparison examples are performed. The impacts of the GPLs volume fraction, distribution type, length-to-depth ratio, humid conditions, external electric voltage, and piezoelectric properties on the critical buckling temperature are studied.

Introduction

Piezoelectric materials have been extensively used in many important industrial models and smart structures.^1–6 This is because of their effective electro-mechanical coupling properties as well as the capability of converting the main energy types, namely mechanical and electrical energies to each other. Consequently, a lot of articles were performed to investigate the behavior of these materials which made dealing with the functionally graded materials (FGMs) more reliable as stated by Wu and his coworkers. ⁷ El Harti et al. ⁸ discussed the impacts of temperature and porosity on the dynamics of a flexible FG porous beam with bonded piezoelectric materials by using the finite element method and the Euler-Bernoulli theory, where the motion equations are derived using Hamilton principle. Mallek et al. ⁹ analyzed the nonlinear dynamic behavior of piezolaminated functionally graded carbon nanotube-reinforced composite (FG-CNTRC) shell by using the first-order shear deformation theory (FSDT). Phuc and Kim Khue ¹⁰ investigated the mechanical behavior of piezoelectric FGM plates resting on two-parameter elastic foundations under thermal conditions depending on the third-order shear deformation theory together with finite element method. Zenkour and Aljadani ¹¹ studied the electro-mechanical buckling behavior of simply supported functionally-graded piezoelectric (FGP) plates subjected to external electric voltage with the help of a quasi-3D refined plate theory. Sobhy ¹² discussed depending on a four variable shell theory the magneto-electro thermal bending of composite doubly-curved shallow shells consisting of two piezoelectromagnetic face sheets considering a variety of boundary conditions. Moreover, Abazid and Sobhy ¹³ investigated the thermal and electro-mechanical size-dependent bending of simply supported functionally graded piezoelectric (FGP) microplates depending on a refined four-variable shear deformation plate theory together with the modified couple stress theory. Furthermore, Abazid ¹⁴ explored the thermomechanical buckling, free vibration, and wave dispersion in smart piezo-electromagnetic nanoplates under humid conditions by using both the nonlocal strain gradient theory (NSGT) and the sinusoidal shear deformation plate theory.

Graphene, a two-dimensional thick layer of carbon atoms joined in a hexagonal form, has a low mass density in addition to superior electromechanical properties. ¹⁵ It has been considered as an efficient reinforcement of the piezoelectric composite structures because it improves their electromechanical characteristics and stiffness.^16,17 Mao et al. ¹⁸ analyzed the small scale influence on the frequencies of graphene nano-platelets reinforced functionally graded piezoelectric composite micro-plate depending on the nonlocal constitutive relation in addition to von-Karman geometric non-linearity, where the nonlinear equations of motion were solved by the differential quadrature (DQ) method. Moreover, Rafiee et al. ¹⁹ investigated the buckling of graphene epoxy nanocomposite beam and deduced that the buckling load increases with a very small increase of the weight fraction of graphene platelets (GPLs) into the epoxy. Sobhy and Alakel Abazid ²⁰ explored the impact of an axial magnetic field on the free vibration and mechanical buckling responses of an FG graphene-reinforced sandwich deep curved nanobeam with viscoelastic core by using the NSGT. Furthermore, Abazid et al. ²¹ discussed the wave propagation in FG porous nanoplate reinforced by uniformly or non-uniformly distributed GPLs and lying on Pasternak elastic foundation depending on both the NSGT and new quasi 3D plate theory, where the motion equations were derived by employing Hamilton’s principle. In another article, Alakel Abazid ²² studied the thermal buckling of porous GPLs-reinforced nanoplate resting on Pasternak foundation and exposed to magnetic field impacts in humid medium. The displacement field was described by employing the modified Reddy’s plate theory and the governing equations were deduced using the NSGT and the principle of virtual displacement. The literature includes a lot of articles that analyze the different behaviors of piezoelectric FG-GPLs like Refs.^23–29

In the current study, the thermal buckling change of piezoelectric functionally graded graphene platelets (PFGGPL) beam which is subjected to external electric voltage in humid environment is discussed. Depending on the Halpin-Tsai model, the effective Young’s modulus of the nanocomposite beam is evaluated. Also according to the rule of the mixture, the Poisson’s ratio, mass density, and piezoelectric properties are computed considering four FGGPL distribution types. To represent the displacement components, a refined two-unknown beam theory is implemented considering both of the shear deformation and the thickness stretching effect. The stability differential equations are deduced by using the principle of virtual work including thermal, moisture, and electric loads. Next, the accuracy of the obtained buckling temperature is validated by comparing with some examples from other researches. Moreover, the influences of the GPLs volume fraction, distribution type, length-to-depth ratio, humid conditions, external electric voltage, and piezoelectric properties on the critical buckling temperature are all investigated.

Problem formulation

In this work, we consider multi-composite layers FG graphene platelets strengthened piezoelectric beam with length $L$ , width $a$ , and total thickness $\tilde{H}$ as illustrated in Figure 1. The present system is assumed to be subjected to an external electric voltage and hygrothermal conditions.

Figure 1.

Configuration of multilayer functionally graded graphene platelets strengthened piezoelectric beam.

Displacement field

A modified refined two-variable beam theory^30,31 that improved by introducing the thickness stretching effects is presented to describe the displacement components that given by

$$\begin{array}{c}\overline{U}(x,z)=u_0(x)-z\frac{\mathrm{d}{w}_1}{\mathrm{d}x}-\varphi (z)\frac{\mathrm{d}{w}_2}{\mathrm{d}x},\hfill \\ \overline{V}(x,z)=0,\overline{W}(x,z)=w_1(x)+w_2(x)+f\prime (z)w_3(x),\hfill \\ \varphi (z)=z-f(z)\hfill \end{array}$$

(1)

where $u_0$ denotes the component of the axial displacement of the mid-plane, the transverse deflection $\overline{W}$ is divided into components $w_1$ , $w_2$ , and $w_3$ . The first two components indicate the bending, shear displacements. While, the function $w_3(x)$ denotes the stretching participation of the displacement that controls the normal strain influence. The function $f(z)$ represents the configure of the shear stress through the thickness of the beam, and it can be written as

$$f(z)=\frac{z}{1+{(z/\tilde{H})}^2}-\frac{5}{8}\frac{z^3}{{\tilde{H}}^2}.$$

(2)

Based on the displacement field (1), one can obtain the components of the strains as:

$$\begin{array}{cc}\multicolumn{1}{c}{{\epsilon }_{\mathrm{xx}}}& =\frac{\partial \overline{U}}{\partial x}=\frac{\mathrm{d}{u}_0}{\mathrm{d}x}-z\frac{{\mathrm{d}}^2{w}_1}{\mathrm{d}{x}^2}-\varphi (z)\frac{{\mathrm{d}}^2{w}_2}{\mathrm{d}{x}^2},\hfill \\ \multicolumn{1}{c}{{\epsilon }_{\mathrm{xz}}}& =\frac{\partial \overline{W}}{\partial x}+\frac{\partial \overline{U}}{\partial z}=(\frac{\mathrm{d}{w}_2}{\mathrm{d}x}+\frac{\mathrm{d}{w}_3}{\mathrm{d}x})f\prime (z)\hfill \\ \multicolumn{1}{c}{{\epsilon }_{\mathrm{zz}}}& =\frac{\partial \overline{W}}{\partial z}=f″(z)w_3(x).\hfill \end{array}$$

(3)

Constitutive relations

With respect to the piezoelasticity theory,^3,4,32 the constitutive relation for the components of the stresses $\sigma$ and the electric displacements $D$ can be given as

$$\begin{array}{c}\multicolumn{1}{c}{{\left\{\begin{array}{c}{\sigma }_{\mathrm{xx}}\hfill \\ \multicolumn{1}{c}{{\sigma }_{\mathrm{zz}}}\\ \multicolumn{1}{c}{{\sigma }_{\mathrm{xz}}}\end{array}\right\}}^i={\left[\begin{array}{ccc}{F}_{11}\hfill & F_{13}\hfill & 0\hfill \\ \multicolumn{1}{c}{F_{13}}& F_{33}\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & F_{55}\hfill \end{array}\right]}^i\left\{\begin{array}{c}{\epsilon }_{\mathrm{xx}}\hfill \\ \multicolumn{1}{c}{{\epsilon }_{\mathrm{zz}}}\\ \multicolumn{1}{c}{{\epsilon }_{\mathrm{xz}}}\end{array}\right\}}\\ \multicolumn{1}{c}{-{\left\{\begin{array}{c}{F}_{11}^{T}T+F_{11}^{C}C\hfill \\ \multicolumn{1}{c}{F_{33}^{T}T+F_{33}^{C}C}\\ \multicolumn{1}{c}{0}\end{array}\right\}}^i-{\left\{\begin{array}{c}{A}_{31}{E}_z\hfill \\ \multicolumn{1}{c}{A_{33}{E}_z}\\ \multicolumn{1}{c}{A_{15}{E}_x}\end{array}\right\}}^i,}\end{array}$$

(4)

$$\begin{array}{c}\multicolumn{1}{c}{{\left\{\begin{array}{c}{D}_x\hfill \\ \multicolumn{1}{c}{D_z}\end{array}\right\}}^i={\left[\begin{array}{ccc}{A}_{15}\hfill & 0\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& A_{31}\hfill & A_{33}\hfill \end{array}\right]}^i{\left\{\begin{array}{c}{\epsilon }_{\mathrm{xz}}\hfill \\ \multicolumn{1}{c}{{\epsilon }_{\mathrm{xx}}-\beta T-\overline{\kappa }C}\\ \multicolumn{1}{c}{{\epsilon }_{\mathrm{zz}}-\beta T-\overline{\kappa }C}\end{array}\right\}}^i}\\ \multicolumn{1}{c}{+{\left\{\begin{array}{c}{\alpha }_{11}{E}_x\hfill \\ \multicolumn{1}{c}{{\alpha }_{33}{E}_z}\end{array}\right\}}^i,}\end{array}$$

(5)

where $(F_{11},F_{13},F_{55})$ , $(F_{11}^T,F_{33}^T)$ , and $(F_{11}^C,F_{33}^C)$ denote, respectively, the elastic coefficients, the thermal coefficients, and moisture coefficients of the FGGPLs composite beam; $T$ and $C$ denote the applied temperature and moisture, respectively; $A_{\mathrm{im}}$ and ${\alpha }_{\mathrm{ii}}$ represent the coefficients of the piezoelectric and the dielectric, respectively. The elastic coefficients, thermal coefficients, and moisture coefficients of the FGGPLs composite beam are written as

$$\begin{array}{c}{F}_{11}^i=F_{33}^i=\frac{E_i}{1-{\left({\nu }_i\right)}^2},F_{13}^i=\frac{{\nu }_i{E}_i}{1-{\left({\nu }_i\right)}^2},F_{55}^i=\frac{E_i}{2(1+{\nu }_i)},\hfill \\ F_{11}^{\mathrm{iT}}=F_{33}^{\mathrm{iT}}=\frac{E_i{\beta }_i}{1-{\nu }_i},F_{11}^{\mathrm{iC}}=F_{33}^{\mathrm{iC}}=\frac{E_i{\overline{\kappa }}_i}{1-{\nu }_i},\hfill \end{array}$$

(6)

The electric field $E_n$ can take the following form ³²

$$E_n=-\nabla \mathrm{\Psi },$$

(7)

where $\mathrm{\Psi }$ represents the electric potential that can be defined as ³²

$$\mathrm{\Psi }(x,z)=-\cos \left(\frac{\pi z}{\tilde{H}}\right)g(x)+\frac{2z{V}_a}{\tilde{H}},$$

(8)

in which $V_a$ is the external applied voltage and $g(x)$ is the electric potential at the mid-plane of the beam. Substituting equation (8) into equation (7) gives

$$\begin{array}{c}\multicolumn{1}{c}{E_x=-\frac{\partial \mathrm{\Psi }}{\partial x}=\cos \left(\frac{\pi z}{\tilde{H}}\right)g\prime (x),}\\ \multicolumn{1}{c}{E_z=-\frac{\partial \mathrm{\Psi }}{\partial z}=-\frac{2V_a}{\tilde{H}}-\frac{\pi }{\tilde{H}}g(x)\sin \left(\frac{\pi z}{\tilde{H}}\right).}\end{array}$$

(9)

Beam construction

Within the framework of Halpin-Tsai model, the effective Young’s moduli $E_i$ for the $i$ th layer of the composite beam are given by ³³

$$E_i=\frac{E_{\mathrm{pz}}}{8}[\frac{3(1+2{\kappa }_1{\delta }_1^{\mathrm{GP}}{V}_i^{\mathrm{GP}})}{1-{\kappa }_1{V}_i^{\mathrm{GP}}}+\frac{5(1+2{\kappa }_2{\delta }_2^{\mathrm{GP}}{V}_i^{\mathrm{GP}})}{1-{\kappa }_2{V}_i^{\mathrm{GP}}}],$$

(10)

where $E_{\mathrm{pz}}$ stands for the piezoelectric Young’s modulus, ${\delta }_1^{\mathrm{GP}}=L^{\mathrm{GP}}/{\tilde{H}}^{\mathrm{GP}}$ and ${\delta }_2^{\mathrm{GP}}=a^{\mathrm{GP}}/{\tilde{H}}^{\mathrm{GP}}$ , in which $L^{\mathrm{GP}}$ , $a^{\mathrm{GP}}$ , and ${\tilde{H}}^{\mathrm{GP}}$ indicate the length, width, and thickness of GPLs, respectively; $V_i^{\mathrm{GP}}$ is the volume fraction of graphene for the $i$ th layer, and the coefficients ${\kappa }_1$ and ${\kappa }_2$ are written as ³³

$${\kappa }_1=\frac{E^{\mathrm{GP}}-E_{\mathrm{pz}}}{E^{\mathrm{GP}}+2{\delta }_1^{\mathrm{GP}}{E}_{\mathrm{pz}}},{\kappa }_2=\frac{E^{\mathrm{GP}}-E_{\mathrm{pz}}}{E^{\mathrm{GP}}+2{\delta }_2^{\mathrm{GP}}{E}_{\mathrm{pz}}},$$

(11)

in which $E^{\mathrm{GP}}$ indicates Young’s modulus of the GPLs. While, the effective material properties of the FGGPLs reinforced composite beam for the $i$ th layer, namely Poisson’s ration ${\nu }_i$ , mass density ${\rho }_i$ , thermal expansion coefficient ${\beta }_i$ , and moisture expansion coefficient $\overline{\kappa }$ can be evaluated as

$$\begin{array}{c}\multicolumn{1}{c}{{\nu }_i=V_i^{\mathrm{GP}}{\nu }^{\mathrm{GP}}+(1-V_i^{\mathrm{GP}}){\nu }_{\mathrm{pz}},}\\ \multicolumn{1}{c}{{\rho }_i=V_i^{\mathrm{GP}}{\rho }^{\mathrm{GP}}+(1-V_i^{\mathrm{GP}}){\rho }_{\mathrm{pz}},}\\ \multicolumn{1}{c}{{\beta }_i=V_i^{\mathrm{GP}}{\beta }^{\mathrm{GP}}+(1-V_i^{\mathrm{GP}}){\beta }_{\mathrm{pz}},}\\ \multicolumn{1}{c}{{\overline{\kappa }}_i=V_i^{\mathrm{GP}}{\overline{\kappa }}^{\mathrm{GP}}+(1-V_i^{\mathrm{GP}}){\overline{\kappa }}_{\mathrm{pz}}}\end{array}$$

(12)

in which ${\nu }^{\mathrm{GP}}({\nu }_{\mathrm{pz}})$ , ${\rho }^{\mathrm{GP}}$ $({\rho }_{\mathrm{pz}})$ , ${\beta }^{\mathrm{GP}}$ $({\beta }_{\mathrm{pz}})$ , and ${\overline{\kappa }}^{\mathrm{GP}}$ $({\overline{\kappa }}_{\mathrm{pz}})$ are Poisson’s ratio of the GPLs (piezoelectric), the mass density of the GPLs (piezoelectric), the thermal expansion coefficient of the GPLs (piezoelectric), and the moisture expansion coefficient of the GPLs (piezoelectric), respectively. Furthermore, the remaining coefficients of piezoelectric and dielectric $(A_{31},A_{15})$ and $({\alpha }_{11},{\alpha }_{33})$ , respectively, can be calculated as follows:

$$\begin{array}{cc}\multicolumn{1}{c}{A_{31}}& =V^{\mathrm{GP}}{A}_{31}^{\mathrm{GP}}+(1-V^{\mathrm{GP}})A_{31,\mathrm{pz}}\hfill \\ \multicolumn{1}{c}{A_{15}}& =V^{\mathrm{GP}}{A}_{15}^{\mathrm{GP}}+(1-V^{\mathrm{GP}})A_{15,\mathrm{pz}}\hfill \\ \multicolumn{1}{c}{{\alpha }_{11}}& =V^{\mathrm{GP}}{\alpha }_{11}^{\mathrm{GP}}+(1-V^{\mathrm{GP}}){\alpha }_{11,\mathrm{pz}}\hfill \\ \multicolumn{1}{c}{{\alpha }_{33}}& =V^{\mathrm{GP}}{\alpha }_{33}^{\mathrm{GP}}+(1-V^{\mathrm{GP}}){\alpha }_{33,\mathrm{pz}},\hfill \end{array}$$

(13)

where $A_{31}^{\mathrm{GP}}$ $(A_{31,\mathrm{pz}})$ , $A_{15}^{\mathrm{GP}}$ $(A_{15,\mathrm{pz}})$ , ${\alpha }_{11}^{\mathrm{GP}}$ $({\alpha }_{11,\mathrm{pz}})$ , and ${\alpha }_{33}^{\mathrm{GP}}$ $({\alpha }_{33,\mathrm{pz}})$ define the piezoelectric coefficients of GPLs (piezoelectric), and the dielectric coefficients of GPLs (piezoelectric), respectively. With respect to a modified piece-wise rule, the volume fraction of GPLs will be diversed across the thickness of the beam layers as presented in Figure 1, and in the present analysis four different patterns of FGGP distribution are studied as follows:

Pattern 1

In this pattern, the graphene platelets can be steadily distributed throughout thickness $\tilde{H}$ . As a result, the volume fraction of the $i$ th layer can be expressed by

$$V_i^{\mathrm{GP}}=V_{\mathrm{m}\mathrm{ax}}^{\mathrm{GP}}=\frac{{\rho }_{\mathrm{pz}}{F}^{\mathrm{GP}}}{{\rho }_{\mathrm{pz}}{F}^{\mathrm{GP}}+{\rho }^{\mathrm{GP}}(1-F^{\mathrm{GP}})},\mathrm{for}\mathrm{SD}$$

(14)

in which $F^{\mathrm{GP}}$ describes the weight fraction of GPLs.

Pattern 2

In this particular pattern, the volume fraction of graphene decreases in a monotonic way from its maximum value $V_{\mathrm{m}\mathrm{ax}}^{\mathrm{GP}}$ at the top-surface to minimum value $V_{\mathrm{m}\mathrm{in}}^{\mathrm{GP}}$ at the bottom-surface of the beam. We can indicate the gradient-distributions of graphene by FG-V. Therefore, the volume fraction for each layer can be written as

$$V_i^{\mathrm{GP}}=V_{\mathrm{m}\mathrm{in}}^{\mathrm{GP}}+(V_{\mathrm{m}\mathrm{ax}}^{\mathrm{GP}}-V_{\mathrm{m}\mathrm{in}}^{\mathrm{GP}}){\left(\frac{i-1}{M-1}\right)}^j\mathrm{for}\mathrm{FG}-\mathrm{V},$$

(15)

in which $M$ represents the number of layers of the beam, and $0\le j\le \infty$ is a power-law index.

Pattern 3

The volume fraction of graphene equals $V_{\mathrm{m}\mathrm{in}}^{\mathrm{GP}}$ at the mid-plane of the sheets. While, it has a maximum value at the upper and lower sheets. Consequently, $V_i^{\mathrm{GP}}$ is expressed as

$$\begin{array}{c}{V}_i^{\mathrm{GP}}=V_{\mathrm{m}\mathrm{in}}^{\mathrm{GP}}+(V_{\mathrm{m}\mathrm{ax}}^{\mathrm{GP}}-V_{\mathrm{m}\mathrm{in}}^{\mathrm{GP}}){\left(\frac{|2i-M-1|-1}{M-2}\right)}^j\mathrm{for}\mathrm{FG}-\mathrm{X}.\hfill \end{array}$$

(16)

Pattern 4

In this state, the volume fraction of graphene equals $V_{\mathrm{m}\mathrm{in}}^{\mathrm{GP}}$ at the top and bottom beam surfaces. While, it has a maximum value at the mid-plane of the beam. Consequently, $V_i^{\mathrm{GP}}$ can be expressed as

$$\begin{array}{c}{V}_i^{\mathrm{GP}}=V_{\mathrm{m}\mathrm{in}}^{\mathrm{GP}}+(V_{\mathrm{m}\mathrm{ax}}^{\mathrm{GP}}-V_{\mathrm{m}\mathrm{in}}^{\mathrm{GP}}){\left(\frac{|2i-M-1|+1-M}{2-M}\right)}^j,\mathrm{for}\mathrm{FG}-\mathrm{O}.\hfill \end{array}$$

(17)

Temperature and moisture fields

The applied temperature $T$ and moisture $C$ are distributed in the thickness direction of the FGGPLs beams, and they can be obtained for each layer depending on the equations of the heat conduction and moisture diffusion as

$$\frac{\mathrm{d}}{\mathrm{d}z}\left({\overline{\mathcal{G}}}_i\frac{\mathrm{d}{\mathrm{\Upsilon }}_i(z)}{\mathrm{d}z}\right)=0,\mathrm{\Upsilon }=T,C,$$

(18)

where ${\overline{\mathcal{G}}}_i$ is the thermal conductivity and moisture diffusivity for the $i$ th layer. Equation (18) will be integrated subject to the boundary conditions of temperature and moisture together with the continuity conditions, which take the forms

$$T(-\frac{\tilde{H}}{2})=T^0,T\left(\frac{\tilde{H}}{2}\right)=T^1,C(-\frac{\tilde{H}}{2})=C^0,C\left(\frac{\tilde{H}}{2}\right)=C^1,{\mathrm{\Upsilon }}_i({\tilde{H}}_{i+1})={\mathrm{\Upsilon }}_{i+1}({\tilde{H}}_{i+1})=T_i^l,\frac{\mathrm{d}{\mathrm{\Upsilon }}_{\mathrm{i}}(z)}{\mathrm{d}z}=\frac{\mathrm{d}{\mathrm{\Upsilon }}_{\mathrm{i}+1}(z)}{\mathrm{d}z}\mathrm{on}z={\tilde{H}}_{i+1},\mathrm{\Upsilon }=T,C,$$

(19)

where ${\tilde{H}}_{i+1}$ indicates the interface coordinate between the layers $i$ and $i+1$ . In addition, the constants $T^0$ $(C^0)$  and $T^1$ $(C^1)$ are the temperature (moisture) of the lower- and upper-surfaces of the beam, respectively. Integration of (18) subject to (19) leads to

$${\mathrm{\Upsilon }}_i(z)=\frac{({\mathrm{\Upsilon }}_i^0-{\mathrm{\Upsilon }}_{i-1}^0)(z-{\tilde{H}}_i)}{{\tilde{H}}_{i+1}-{\tilde{H}}_i}+{\mathrm{\Upsilon }}_{i-1}^0,\mathrm{\Upsilon }=T,C,$$

(20)

where $(i=1,2,\dots p)$ denotes the number of the structure layers, ${\mathrm{\Upsilon }}_0^0={\mathrm{\Upsilon }}^0$ , ${\mathrm{\Upsilon }}_p^1={\mathrm{\Upsilon }}^1$ . Equation (20) gives the general mechanism of the temperature and moisture.

Governing equations

The equations of motion associated with the displacement field in equation (1) can be derived by utilizing the variation of the strain energy $\delta U_{\mathrm{S}\mathrm{E}}$ and the work done by the external force $\delta U_{\mathrm{F}\mathrm{E}}$ which can be stated by

$$\begin{array}{c}\begin{array}{c}\delta U_{\mathrm{S}E}=\sum _{i=1}^M{\int }_{{\tilde{H}}_i}^{{\tilde{H}}_{i+1}}{\int }_0^L({\sigma }_{\mathrm{xx}}^i\delta {\epsilon }_{\mathrm{xx}}+{\sigma }_{\mathrm{xz}}^i\delta {\epsilon }_{\mathrm{xz}}+{\sigma }_{\mathrm{zz}}^i\delta {\epsilon }_{\mathrm{zz}}\hfill \\ \multicolumn{1}{c}{-D_x^i\delta E_x-D_z^i\delta E_z)\mathrm{d}x\mathrm{d}z,}\end{array}\hfill \\ \begin{array}{c}\delta U_{\mathrm{F}E}={\int }_0^L(G_1^T+G_1^C+G_1^E)(\frac{{\mathrm{d}}^2{w}_1}{\mathrm{d}{x}^2}+\frac{{\mathrm{d}}^2{w}_2}{\mathrm{d}{x}^2}+\frac{{\mathrm{d}}^2{w}_3}{\mathrm{d}{x}^2})\hfill \\ \multicolumn{1}{c}{\delta (w_1+w_2+w_3)\mathrm{d}x,}\end{array}\hfill \end{array}$$

(21)

in which $G_1^T$ , $G_1^C$ , and $G_1^E$ are the in-plane hygrothermal forces which are given, respectively, as

$$\begin{array}{c}{G}_1^T=-\sum _{i=1}^M{\int }_{{\tilde{H}}_i}^{{\tilde{H}}_{i+1}}{F}_{11}^{\mathrm{iT}}T\mathrm{d}z,G_1^C=-\sum _{i=1}^M{\int }_{{\tilde{H}}_i}^{{\tilde{H}}_{i+1}}{F}_{11}^{\mathrm{iC}}C\mathrm{d}z,\hfill \\ G_1^E=\sum _{i=1}^M{\int }_{{\tilde{H}}_i}^{{\tilde{H}}_{i+1}}{A}_{31}^i\left(\frac{2V_a}{\tilde{H}}\right)\mathrm{d}z,\hfill \end{array}$$

(22)

For more extension of equation (21), equation (3) is substituted into equation (21) to yield

$$\begin{array}{l}\delta U_{\mathrm{SE}}={\int }_0^L({\mathcal{N}}_x\delta \frac{\mathrm{d}{u}_0}{\mathrm{d}x}-{ℳ}_x\delta \frac{{\mathrm{d}}^2{w}_1}{\mathrm{d}{x}^2}-{ℱ}_x\delta \frac{{\mathrm{d}}^2{w}_2}{\mathrm{d}{x}^2}\\ +{\mathcal{Q}}_x\delta \left(\frac{\mathrm{d}{w}_2}{\mathrm{d}x}+\frac{\mathrm{d}{w}_3}{\mathrm{d}x}\right)\\ +{\mathscr{H}}_x\delta w_3-{ℛ}_x\delta \frac{\partial g}{\partial x}+{ℛ}_z\delta g)\mathrm{d}x,\end{array}$$

(23)

where

$$\begin{array}{c}{\mathcal{N}}_x=\sum _{i=1}^M{\int }_{{\tilde{H}}_i}^{{\tilde{H}}_{i+1}}{\sigma }_{\mathrm{xx}}^i\mathrm{d}z,{\mathcal{M}}_x=\sum _{i=1}^M{\int }_{{\tilde{H}}_i}^{{\tilde{H}}_{i+1}}z{\sigma }_{\mathrm{xx}}^i\mathrm{d}z,\hfill \\ {\mathcal{F}}_x=\sum _{i=1}^M{\int }_{{\tilde{H}}_i}^{{\tilde{H}}_{i+1}}\varphi (z){\sigma }_{\mathrm{xx}}^i\mathrm{d}z,{\mathcal{Q}}_x=\sum _{i=1}^M{\int }_{{\tilde{H}}_i}^{{\tilde{H}}_{i+1}}f\prime (z){\sigma }_{\mathrm{xz}}^i\mathrm{d}z,\hfill \\ {\mathcal{H}}_x=\sum _{i=1}^M{\int }_{{\tilde{H}}_i}^{{\tilde{H}}_{i+1}}f″(z){\sigma }_{\mathrm{zz}}^i\mathrm{d}z,\hfill \\ {\mathcal{R}}_x=\sum _{i=1}^M{\int }_{{\tilde{H}}_i}^{{\tilde{H}}_{i+1}}{D}_x^i\cos \left(\frac{\pi z}{\tilde{H}}\right)\mathrm{d}z,\hfill \\ {\mathcal{R}}_z=\sum _{i=1}^M{\int }_{{\tilde{H}}_i}^{{\tilde{H}}_{i+1}}{D}_z^i\frac{\pi }{\tilde{H}}\sin \left(\frac{\pi z}{\tilde{H}}\right)\mathrm{d}z.\hfill \end{array}$$

(24)

The stability equations can be deduced from the principle of virtual work which is:

$$\begin{array}{c}\delta U_{\mathrm{S}\mathrm{E}}-\delta U_{\mathrm{F}\mathrm{E}}=0\hfill \end{array}$$

(25)

Substituting equations (21) and (23) into equation (25) gives the stability equations as:

$$\begin{array}{c}\frac{\mathrm{d}{\mathcal{N}}_x}{\mathrm{d}x}=0,\hfill \\ \frac{{\mathrm{d}}^2{\mathcal{M}}_x}{\mathrm{d}{x}^2}+(G_1^T+G_1^C+G_1^E)\hfill \\ (\frac{{\mathrm{d}}^2{w}_1}{\mathrm{d}{x}^2}+\frac{{\mathrm{d}}^2{w}_2}{\mathrm{d}{x}^2}+\frac{{\mathrm{d}}^2{w}_3}{\mathrm{d}{x}^2})=0,\hfill \\ \frac{{\mathrm{d}}^2{\mathcal{F}}_x}{\mathrm{d}{x}^2}+\frac{\mathrm{d}{\mathcal{Q}}_x}{\mathrm{d}x}+(G_1^T+G_1^C+G_1^E)\hfill \\ (\frac{{\mathrm{d}}^2{w}_1}{\mathrm{d}{x}^2}+\frac{{\mathrm{d}}^2{w}_2}{\mathrm{d}{x}^2}+\frac{{\mathrm{d}}^2{w}_3}{\mathrm{d}{x}^2})=0,\hfill \\ \frac{\mathrm{d}{\mathcal{Q}}_x}{\mathrm{d}x}-{\mathcal{H}}_x+(G_1^T+G_1^C+G_1^E)\hfill \\ (\frac{{\mathrm{d}}^2{w}_1}{\mathrm{d}{x}^2}+\frac{{\mathrm{d}}^2{w}_2}{\mathrm{d}{x}^2}+\frac{{\mathrm{d}}^2{w}_3}{\mathrm{d}{x}^2})=0,\hfill \\ \frac{\mathrm{d}{\mathcal{R}}_x}{\mathrm{d}x}+{\mathcal{R}}_z=0.\hfill \end{array}$$

(26)

By inserting equations (4) and (5) into equation (24) with the help of equations (3) and (9), one obtains the stress resultants as:

$$\begin{array}{cc}\multicolumn{1}{c}{{\mathcal{N}}_x=}& K_{11}\frac{\mathrm{d}{u}_0}{\mathrm{d}x}-K_{12}\frac{{\mathrm{d}}^2{w}_1}{\mathrm{d}{x}^2}-K_{13}\frac{{\mathrm{d}}^2{w}_2}{\mathrm{d}{x}^2}+K_{14}{w}_3+K_{15}g+G_1^T+G_1^C+G_1^E,\hfill \\ \multicolumn{1}{c}{{\mathcal{M}}_x=}& K_{12}\frac{\mathrm{d}{u}_0}{\mathrm{d}x}-K_{22}\frac{{\mathrm{d}}^2{w}_1}{\mathrm{d}{x}^2}-K_{23}\frac{{\mathrm{d}}^2{w}_2}{\mathrm{d}{x}^2}+K_{24}{w}_3+K_{25}g+{\mathcal{M}}_x^T+{\mathcal{M}}_x^C+{\mathcal{M}}_x^E,\hfill \\ \multicolumn{1}{c}{{\mathcal{F}}_x=}& K_{13}\frac{\mathrm{d}{u}_0}{\mathrm{d}x}-K_{23}\frac{{\mathrm{d}}^2{w}_1}{\mathrm{d}{x}^2}-K_{33}\frac{{\mathrm{d}}^2{w}_2}{\mathrm{d}{x}^2}+K_{34}{w}_3+K_{35}g+{\mathcal{F}}_x^T+{\mathcal{F}}_x^C+{\mathcal{F}}_x^E,\hfill \\ \multicolumn{1}{c}{{\mathcal{Q}}_x=}& K_{41}(\frac{\mathrm{d}{w}_2}{\mathrm{d}x}+\frac{\mathrm{d}{w}_3}{\mathrm{d}x})-K_{42}g\prime (x),\hfill \\ \multicolumn{1}{c}{{\mathcal{H}}_x=}& K_{14}\frac{\mathrm{d}{u}_0}{\mathrm{d}x}-K_{24}\frac{{\mathrm{d}}^2{w}_1}{\mathrm{d}{x}^2}-K_{34}\frac{{\mathrm{d}}^2{w}_2}{\mathrm{d}{x}^2}+K_{54}{w}_3+K_{55}g+{\mathcal{H}}_x^T+{\mathcal{H}}_x^C+{\mathcal{H}}_x^E,\hfill \\ \multicolumn{1}{c}{{\mathcal{R}}_x=}& K_{42}(\frac{\mathrm{d}{w}_2}{\mathrm{d}x}+\frac{\mathrm{d}{w}_3}{\mathrm{d}x})-K_{62}g\prime (x),\hfill \\ \multicolumn{1}{c}{{\mathcal{R}}_z=}& K_{15}\frac{\mathrm{d}{u}_0}{\mathrm{d}x}-K_{25}\frac{{\mathrm{d}}^2{w}_1}{\mathrm{d}{x}^2}-K_{35}\frac{{\mathrm{d}}^2{w}_2}{\mathrm{d}{x}^2}+K_{55}{w}_3+K_{75}g+{\mathcal{R}}_z^T+{\mathcal{R}}_z^C+{\mathcal{H}}_z^E,\hfill \end{array}$$

(27)

where

$$\begin{array}{c}{K}_{11}=\sum _{i=1}^M{\int }_{{\tilde{H}}_i}^{{\tilde{H}}_{i+1}}{F}_{11}^i\mathrm{d}z,K_{12}=\sum _{i=1}^M{\int }_{{\tilde{H}}_i}^{{\tilde{H}}_{i+1}}{\mathrm{zF}}_{11}^i\mathrm{d}z,K_{13}=\sum _{i=1}^M{\int }_{{\tilde{H}}_i}^{{\tilde{H}}_{i+1}}\varphi (z)F_{11}^i\mathrm{d}z\hfill \\ K_{14}=\sum _{i=1}^M{\int }_{{\tilde{H}}_i}^{{\tilde{H}}_{i+1}}f″F_{13}^i\mathrm{d}z,K_{15}=\sum _{i=1}^M{\int }_{{\tilde{H}}_i}^{{\tilde{H}}_{i+1}}{A}_{31}^i\left(\frac{\pi }{\tilde{H}}\right)\sin \left(\frac{\pi z}{\tilde{H}}\right)\mathrm{d}z,\hfill \\ K_{22}=\sum _{i=1}^M{\int }_{{\tilde{H}}_i}^{{\tilde{H}}_{i+1}}{z}^2{F}_{11}^i\mathrm{d}z,K_{23}=\sum _{i=1}^M{\int }_{{\tilde{H}}_i}^{{\tilde{H}}_{i+1}}z\varphi (z)F_{11}^i\mathrm{d}z\hfill \\ K_{24}=\sum _{i=1}^M{\int }_{{\tilde{H}}_i}^{{\tilde{H}}_{i+1}}\mathrm{zf}″F_{13}^i\mathrm{d}z,K_{25}=\sum _{i=1}^M{\int }_{{\tilde{H}}_i}^{{\tilde{H}}_{i+1}}{\mathrm{zA}}_{31}^i\left(\frac{\pi }{\tilde{H}}\right)\sin \left(\frac{\pi z}{\tilde{H}}\right)\mathrm{d}z,\hfill \\ {\mathcal{M}}_x^E=\sum _{i=1}^M{\int }_{{\tilde{H}}_i}^{{\tilde{H}}_{i+1}}{\mathrm{zA}}_{31}^i\left(\frac{2V_a}{\tilde{H}}\right)\mathrm{d}z,{\mathcal{M}}_x^T=-\sum _{i=1}^M{\int }_{{\tilde{H}}_i}^{{\tilde{H}}_{i+1}}{\mathrm{zF}}_{11}^{\mathrm{iT}}T\mathrm{d}z,\hfill \\ {\mathcal{M}}_x^C=-\sum _{i=1}^M{\int }_{{\tilde{H}}_i}^{{\tilde{H}}_{i+1}}{\mathrm{zF}}_{11}^{\mathrm{iC}}C\mathrm{d}z,K_{33}=\sum _{i=1}^M{\int }_{{\tilde{H}}_i}^{{\tilde{H}}_{i+1}}\varphi {(z)}^2{F}_{11}^i\mathrm{d}z\hfill \\ K_{34}=\sum _{i=1}^M{\int }_{{\tilde{H}}_i}^{{\tilde{H}}_{i+1}}f″\varphi (z)F_{13}^i\mathrm{d}z,K_{35}=\sum _{i=1}^M{\int }_{{\tilde{H}}_i}^{{\tilde{H}}_{i+1}}{A}_{31}^i\varphi (z)\left(\frac{\pi }{\tilde{H}}\right)\sin \left(\frac{\pi z}{\tilde{H}}\right)\mathrm{d}z,\hfill \\ {\mathcal{F}}_x^E=\sum _{i=1}^M{\int }_{{\tilde{H}}_i}^{{\tilde{H}}_{i+1}}{A}_{31}^i\varphi (z)\left(\frac{2V_a}{\tilde{H}}\right)\mathrm{d}z,{\mathcal{F}}_x^T=-\sum _{i=1}^M{\int }_{{\tilde{H}}_i}^{{\tilde{H}}_{i+1}}\varphi (z)F_{11}^{\mathrm{iT}}T\mathrm{d}z,\hfill \\ {\mathcal{F}}_x^C=-\sum _{i=1}^M{\int }_{{\tilde{H}}_i}^{{\tilde{H}}_{i+1}}\varphi (z)F_{11}^{\mathrm{iC}}C\mathrm{d}z,K_{41}=\sum _{i=1}^M{\int }_{{\tilde{H}}_i}^{{\tilde{H}}_{i+1}}{f\prime }^2{F}_{55}^i\mathrm{d}z,\hfill \\ K_{42}=\sum _{i=1}^M{\int }_{{\tilde{H}}_i}^{{\tilde{H}}_{i+1}}{A}_{15}^{i}f\prime \cos \left(\frac{\pi z}{\tilde{H}}\right)\mathrm{d}z,K_{54}=\sum _{i=1}^M{\int }_{{\tilde{H}}_i}^{{\tilde{H}}_{i+1}}{f″}^2{F}_{33}^i\mathrm{d}z,\hfill \\ K_{55}=\sum _{i=1}^M{\int }_{{\tilde{H}}_i}^{{\tilde{H}}_{i+1}}{A}_{33}^{i}f″\left(\frac{\pi }{\tilde{H}}\right)\sin \left(\frac{\pi z}{\tilde{H}}\right)\mathrm{d}z,{\mathcal{H}}_x^E=\sum _{i=1}^M{\int }_{{\tilde{H}}_i}^{{\tilde{H}}_{i+1}}{A}_{33}^{i}f″\left(\frac{2V_a}{\tilde{H}}\right)\mathrm{d}z,\hfill \\ {\mathcal{H}}_x^T=-\sum _{i=1}^M{\int }_{{\tilde{H}}_i}^{{\tilde{H}}_{i+1}}f″F_{33}^{\mathrm{iT}}T\mathrm{d}z,{\mathcal{H}}_x^C=-\sum _{i=1}^M{\int }_{{\tilde{H}}_i}^{{\tilde{H}}_{i+1}}f″F_{33}^{\mathrm{iC}}C\mathrm{d}z,\hfill \\ K_{62}=\sum _{i=1}^M{\int }_{{\tilde{H}}_i}^{{\tilde{H}}_{i+1}}{\alpha }_{11}^i\stackrel{2}{\cos }\left(\frac{\pi z}{\tilde{H}}\right)\mathrm{d}z,K_{75}=-\sum _{i=1}^M{\int }_{{\tilde{H}}_i}^{{\tilde{H}}_{i+1}}{\alpha }_{33}^i{\left(\frac{\pi }{\tilde{H}}\right)}^2\stackrel{2}{\sin }\left(\frac{\pi z}{\tilde{H}}\right)\mathrm{d}z,\hfill \\ {\mathcal{R}}_z^E=-\sum _{i=1}^M{\int }_{{\tilde{H}}_i}^{{\tilde{H}}_{i+1}}{\alpha }_{33}^i\left(\frac{2V_a}{\tilde{H}}\right)\left(\frac{\pi }{\tilde{H}}\right)\sin \left(\frac{\pi z}{\tilde{H}}\right)\mathrm{d}z,\hfill \\ {\mathcal{R}}_z^T=-\sum _{i=1}^M{\int }_{{\tilde{H}}_i}^{{\tilde{H}}_{i+1}}(A_{33}^i+A_{31}^i){\beta }_{i}T\left(\frac{\pi }{\tilde{H}}\right)\sin \left(\frac{\pi z}{\tilde{H}}\right)\mathrm{d}z,\hfill \\ {\mathcal{R}}_z^C=-\sum _{i=1}^M{\int }_{{\tilde{H}}_i}^{{\tilde{H}}_{i+1}}(A_{33}^i+A_{31}^i){\overline{\kappa }}_{i}C\left(\frac{\pi }{\tilde{H}}\right)\sin \left(\frac{\pi z}{\tilde{H}}\right)\mathrm{d}z.\hfill \end{array}$$

(28)

Substituting equation (27) into equation (26) leads to the governing equations in terms of the displacement components ( $u_0,w_1,w_2,w_3$ ) and the electric potential ( $g$ ) as:

$$\begin{array}{c}{K}_{11}\frac{{\mathrm{d}}^2{u}_0}{\mathrm{d}{x}^2}-K_{12}\frac{{\mathrm{d}}^3{w}_1}{\mathrm{d}{x}^3}-K_{13}\frac{{\mathrm{d}}^3{w}_2}{\mathrm{d}{x}^3}+K_{14}\frac{\mathrm{d}{w}_3}{\mathrm{d}x}+K_{15}\frac{\mathrm{d}g}{\mathrm{d}x}=0,\hfill \\ \mathrm{\Delta }{\mathrm{TK}}_{12}^T\frac{{\mathrm{d}}^2{w}_3}{\mathrm{d}{x}^2}+\mathrm{\Delta }{\mathrm{TK}}_{12}^T\frac{{\mathrm{d}}^2{w}_1}{\mathrm{d}{x}^2}+\mathrm{\Delta }{\mathrm{TK}}_{12}^T\frac{{\mathrm{d}}^2{w}_2}{\mathrm{d}{x}^2}+T^0{K}_{11}^T\frac{{\mathrm{d}}^2{w}_3}{\mathrm{d}{x}^2}+T^0{K}_{11}^T\frac{{\mathrm{d}}^2{w}_1}{\mathrm{d}{x}^2}+T^0{K}_{11}^T\frac{{\mathrm{d}}^2{w}_2}{\mathrm{d}{x}^2}+G_1^C\frac{{\mathrm{d}}^2{w}_3}{\mathrm{d}{x}^2}+\hfill \\ G_1^C\frac{{\mathrm{d}}^2{w}_1}{\mathrm{d}{x}^2}+G_1^C\frac{{\mathrm{d}}^2{w}_2}{\mathrm{d}{x}^2}+G_1^E\frac{{\mathrm{d}}^2{w}_3}{\mathrm{d}{x}^2}+G_1^E\frac{{\mathrm{d}}^2{w}_1}{\mathrm{d}{x}^2}+G_1^E\frac{{\mathrm{d}}^2{w}_2}{\mathrm{d}{x}^2}+K_{12}\frac{{\mathrm{d}}^3{u}_0}{\mathrm{d}{x}^3}-K_{22}\frac{{\mathrm{d}}^4{w}_1}{\mathrm{d}{x}^4}-K_{23}\frac{{\mathrm{d}}^4{w}_2}{\mathrm{d}{x}^4}+\hfill \\ K_{24}\frac{{\mathrm{d}}^2{w}_3}{\mathrm{d}{x}^2}+K_{25}\frac{{\mathrm{d}}^{2}g}{\mathrm{d}{x}^2}=0,\hfill \\ \mathrm{\Delta }T\frac{{\mathrm{d}}^2{w}_3}{\mathrm{d}{x}^2}{K}_{12}^T+\mathrm{\Delta }T\frac{{\mathrm{d}}^2{w}_1}{\mathrm{d}{x}^2}{K}_{12}^T+\mathrm{\Delta }T\frac{{\mathrm{d}}^2{w}_2}{\mathrm{d}{x}^2}{K}_{12}^T+\frac{{\mathrm{d}}^2{w}_3}{\mathrm{d}{x}^2}{T}^0{K}_{11}^T+\frac{{\mathrm{d}}^2{w}_1}{\mathrm{d}{x}^2}{T}^0{K}_{11}^T+\frac{{\mathrm{d}}^2{w}_2}{\mathrm{d}{x}^2}{T}^0{K}_{11}^T+G_1^C\frac{{\mathrm{d}}^2{w}_3}{\mathrm{d}{x}^2}+\hfill \\ G_1^C\frac{{\mathrm{d}}^2{w}_1}{\mathrm{d}{x}^2}+G_1^C\frac{{\mathrm{d}}^2{w}_2}{\mathrm{d}{x}^2}+G_1^E\frac{{\mathrm{d}}^2{w}_3}{\mathrm{d}{x}^2}+G_1^E\frac{{\mathrm{d}}^2{w}_1}{\mathrm{d}{x}^2}+G_1^E\frac{{\mathrm{d}}^2{w}_2}{\mathrm{d}{x}^2}+K_{13}\frac{{\mathrm{d}}^3{u}_0}{\mathrm{d}{x}^3}-K_{23}\frac{{\mathrm{d}}^4{w}_1}{\mathrm{d}{x}^4}-K_{33}\frac{{\mathrm{d}}^4{w}_2}{\mathrm{d}{x}^4}+\hfill \\ K_{34}\frac{{\mathrm{d}}^2{w}_3}{\mathrm{d}{x}^2}+K_{41}\frac{{\mathrm{d}}^2{w}_3}{\mathrm{d}{x}^2}+K_{35}\frac{{\mathrm{d}}^{2}g}{\mathrm{d}{x}^2}-K_{42}\frac{{\mathrm{d}}^{2}g}{\mathrm{d}{x}^2}+K_{41}\frac{{\mathrm{d}}^2{w}_2}{\mathrm{d}{x}^2}=0,\hfill \\ \mathrm{\Delta }T\frac{{\mathrm{d}}^2{w}_3}{\mathrm{d}{x}^2}{K}_{12}^T+\mathrm{\Delta }T\frac{{\mathrm{d}}^2{w}_1}{\mathrm{d}{x}^2}{K}_{12}^T+\mathrm{\Delta }T\frac{{\mathrm{d}}^2{w}_2}{\mathrm{d}{x}^2}{K}_{12}^T+\frac{{\mathrm{d}}^2{w}_3}{\mathrm{d}{x}^2}{T}^0{K}_{11}^T+\frac{{\mathrm{d}}^2{w}_1}{\mathrm{d}{x}^2}{T}^0{K}_{11}^T+\frac{{\mathrm{d}}^2{w}_2}{\mathrm{d}{x}^2}{T}^0{K}_{11}^T+G_1^C\frac{{\mathrm{d}}^2{w}_3}{\mathrm{d}{x}^2}+\hfill \\ G_1^C\frac{{\mathrm{d}}^2{w}_1}{\mathrm{d}{x}^2}+G_1^C\frac{{\mathrm{d}}^2{w}_2}{\mathrm{d}{x}^2}+G_1^E\frac{{\mathrm{d}}^2{w}_3}{\mathrm{d}{x}^2}+G_1^E\frac{{\mathrm{d}}^2{w}_1}{\mathrm{d}{x}^2}+G_1^E\frac{{\mathrm{d}}^2{w}_2}{\mathrm{d}{x}^2}+K_{41}\frac{{\mathrm{d}}^2{w}_3}{\mathrm{d}{x}^2}-K_{42}\frac{{\mathrm{d}}^{2}g}{\mathrm{d}{x}^2}+K_{24}\frac{{\mathrm{d}}^2{w}_1}{\mathrm{d}{x}^2}+\hfill \\ K_{34}\frac{{\mathrm{d}}^2{w}_2}{\mathrm{d}{x}^2}+K_{41}\frac{{\mathrm{d}}^2{w}_2}{\mathrm{d}{x}^2}-K_{14}\frac{\mathrm{d}{u}_0}{\mathrm{d}x}-K_{54}{w}_3(x)-K_{55}g(x)=0,\hfill \\ -K_{42}\frac{{\mathrm{d}}^2{w}_3}{\mathrm{d}{x}^2}+K_{62}\frac{{\mathrm{d}}^{2}g}{\mathrm{d}{x}^2}+K_{25}\frac{{\mathrm{d}}^2{w}_1}{\mathrm{d}{x}^2}+K_{35}\frac{{\mathrm{d}}^2{w}_2}{\mathrm{d}{x}^2}-K_{42}\frac{{\mathrm{d}}^2{w}_2}{\mathrm{d}{x}^2}-K_{15}\frac{\mathrm{d}{u}_0}{\mathrm{d}x}-K_{55}{w}_3+K_{75}g=0,\hfill \end{array}$$

(29)

where

$$\begin{array}{c}{K}_{11}^T=-\sum _{i=1}^M{\int }_{{\tilde{H}}_i}^{{\tilde{H}}_{i+1}}{F}_{11}^{\mathrm{iT}}\mathrm{d}z,\hfill \\ K_{12}^T=-\sum _{i=1}^M{\int }_{{\tilde{H}}_i}^{{\tilde{H}}_{i+1}}{F}_{11}^{\mathrm{iT}}\frac{(z-{\tilde{H}}_i)}{{\tilde{H}}_{i+1}-{\tilde{H}}_i}\mathrm{d}z,\hfill \\ \mathrm{\Delta }T=T^1-T^0.\hfill \end{array}$$

(30)

Thermal buckling solution

The analytical solution of equation (29) is presented according to Navier’s solution procedure. The simply supported boundary conditions for the current shear and normal deformations beam theory at $x=0$ and $x=L$ are given as:

$$\begin{array}{c}{w}_1=w_2=w_3={\mathcal{N}}_x={\mathcal{M}}_x={\mathcal{F}}_x={\mathcal{H}}_x=0.\hfill \end{array}$$

(31)

Accordingly, the displacements that satisfy the above boundary conditions can be given as:

$$\begin{array}{c}{u}_0=U_0\cos (\lambda x),\hfill \\ \{w_1,w_2,w_3,g\}=\{W_1,W_2,W_3,G_0\}\sin (\lambda x),\hfill \end{array}$$

(32)

where $\lambda =n\pi /L$ , $n$ is the mode number, $U_0,W_1,W_2,W_3$ , and $G_0$ are unknown parameters that can be determined. Inserting equation (32) into equation (29) gives the stability equations in the following form:

$$\begin{array}{c}\left[\mathcal{P}\right]\{\mathrm{\Psi }\}=0,\hfill \end{array}$$

(33)

where $\mathrm{\Psi }$ denotes the columns

$$\begin{array}{c}\{\mathrm{\Psi }\}=\{U_0{W}_1{W}_2{W}_3{G}_0\},\hfill \end{array}$$

(34)

and the elements ${\mathcal{P}}_{\mathrm{ij}}={\mathcal{P}}_{\mathrm{ji}}$ of matrix $[\mathcal{P}]$ are expressed as:

$$\begin{array}{c}{\mathcal{P}}_{11}=-{\lambda }^2{K}_{11},{\mathcal{P}}_{12}={\lambda }^3{K}_{12},{\mathcal{P}}_{13}={\lambda }^3{K}_{13},\hfill \\ {\mathcal{P}}_{14}=\lambda K_{14},{\mathcal{P}}_{15}=\lambda K_{15},\hfill \\ {\mathcal{P}}_{22}=-{\lambda }^2({\lambda }^2{K}_{22}+\mathrm{\Delta }{\mathrm{TK}}_{12}^T+T^0{K}_{11}^T+G_1^C+G_1^E),\hfill \\ {\mathcal{P}}_{23}=-{\lambda }^2({\lambda }^2{K}_{23}+\mathrm{\Delta }{\mathrm{TK}}_{12}^T+T^0{K}_{11}^T+G_1^C+G_1^E),\hfill \\ {\mathcal{P}}_{24}=-{\lambda }^2(K_{24}+\mathrm{\Delta }{\mathrm{TK}}_{12}^T+T^0{K}_{11}^T+G_1^C+G_1^E),\hfill \\ {\mathcal{P}}_{25}=-{\lambda }^2{K}_{25},\hfill \\ {\mathcal{P}}_{33}=-{\lambda }^2({\lambda }^2{K}_{33}+K_{41}+\mathrm{\Delta }{\mathrm{TK}}_{12}^T+T^0{K}_{11}^T+G_1^C+G_1^E),\hfill \\ {\mathcal{P}}_{34}=-{\lambda }^2(K_{34}+K_{41}+\mathrm{\Delta }{\mathrm{TK}}_{12}^T+T^0{K}_{11}^T+G_1^C+G_1^E),\hfill \\ {\mathcal{P}}_{35}=-{\lambda }^2(K_{35}-K_{42}),\hfill \\ {\mathcal{P}}_{44}=-{\lambda }^2(K_{41}+\mathrm{\Delta }{\mathrm{TK}}_{12}^T+T^0{K}_{11}^T+G_1^C+G_1^E)-K_{54},\hfill \\ {\mathcal{P}}_{45}={\lambda }^2{K}_{42}-K_{55},\hfill \\ {\mathcal{P}}_{55}={\lambda }^2{K}_{62}+K_{75}.\hfill \end{array}$$

(35)

For non-trivial solution of equation (33), the determinant $|\mathcal{P}|$ should be equal zero. By solving the equation $|\mathcal{P}|=0$ , one can obtain the buckling temperature change $\mathrm{\Delta }T$ .

Numerical results

Numerical results for critical buckling temperature of FGGPL strengthened piezoelectric beams subjected to external electric voltage and humid conditions are presented in this section. The beam thickness is taken as $\tilde{H}=0.004$ m. In addition, the following data are used in the current analysis (unless otherwise stated): $F^{GP}=0.1,e_0=300\times {10}^3,\mathrm{\Delta }C=0.01\%,\gamma =10,E_0=0.01$ , $L^{\mathrm{GP}}=15$  nm, $a^{\mathrm{GP}}=9$  nm, ${\tilde{H}}^{\mathrm{GP}}=0.188$  nm, $T^0=100$  K, $C^0=0$ %. The properties of the piezoelectric material are defined as ² : $E_{\mathrm{pz}}=1.4$  GPa, ${\nu }_{\mathrm{pz}}=0.29,{\rho }_{\mathrm{pz}}=1920$  kg/m³, $A_{31,\mathrm{pz}}=50.535\times {10}^{-3}$  C/m², $A_{33,\mathrm{pz}}=13.212\times {10}^{-3}$  C/m², $A_{15,\mathrm{pz}}=-15.93\times {10}^{-3}$  C/m², ${\alpha }_{11,\mathrm{pz}}=0.5385\times {10}^{-9}$  C/Vm, ${\alpha }_{33,\mathrm{pz}}=0.59571\times {10}^{-9}$  C/Vm, ${\beta }_{\mathrm{pz}}=60\times {10}^{-6}$  K, ${\overline{\kappa }}_{\mathrm{pz}}=0.44$ ${(\mathrm{wt}.\%{\mathrm{H}}_2\mathrm{O})}^{-1}$ . Whereas, the GPL properties are defined as ² : $E^{\mathrm{GP}}=1.01$  TPa, ${\rho }^{\mathrm{GP}}=1060$  kg/m³, ${\nu }^{\mathrm{GP}}=0.186$ , $A_{31}^{\mathrm{GP}}=e_0{A}_{31,\mathrm{pz}}$ , $A_{33}^{\mathrm{GP}}=e_0{A}_{33,\mathrm{pz}}$ , $A_{15}^{\mathrm{GP}}=e_0{A}_{15,\mathrm{pz}}$ , ${\alpha }_{11}^{\mathrm{GP}}=e_0{\alpha }_{11,\mathrm{pz}}$ , ${\alpha }_{33}^{\mathrm{GP}}=e_0{\alpha }_{33,\mathrm{pz}}$ , ${\beta }^{\mathrm{GP}}=5\times {10}^{-6}$  K, ${\overline{\kappa }}^{\mathrm{GP}}=0.26\times {10}^{-3}$ ${(\mathrm{wt}.\%{\mathrm{H}}_2\mathrm{O})}^{-1}$ , where $e_0$ is named as the piezoelectric multiple. ²

A comparison between the obtained critical buckling temperature of functionally graded beam and that obtained by Kiani and Eslami ³⁴ is investigated in Figure 2. The results of Kiani and Eslami ³⁴ has been obtained by employing the Euler-Bernoulli beam theory (EBT) that neglects the shear deformation effect. While, the refined two-variable shear deformation beam theory (RSDT) is utilized to calculate the present results. As expected for the thick beam, a noticeable difference between the results of the EBT ³⁴ (the classical beam theory) and those obtained by the RSDT is obvious, indicating the shear deformation effects. Indeed, the effect of the shear strain is small for thin beams; while, it increases as the beams thickness increase. Therefore, the present theory including the shear deformation, predicts more accurate results than the EBT, especially for thick beams. Whereas, the results of the two theories are in agreement for a thin beam.

Figure 2.

Comparing the obtained critical buckling temperature with that of Kiani and Eslami ³⁴ for an functionally graded beam ( $\overline{\gamma }=\tilde{H}/L$ ).

Figure 3 shows the changes of the critical buckling temperature $T_{\mathrm{cr}}={10}^{-3}\mathrm{\Delta }T$ of FGGPL strengthened piezoelectric beams with respect to various values of external electric voltage $E_0$ versus the length-to-depth ratio $\gamma =L/\tilde{H}$ . Here, the four patterns 1, 2, 3, and 4 of FGGP are considered in the corresponding figures a, b, c, and d, respectively. It is noted that increasing the electric voltage enhances the strength of the nanocomposite beams. Therefore, the critical buckling temperature $T_{\mathrm{cr}}$ increases by increasing the external electric voltage $E_0$ . While, with decreasing the beam thickness, the beam becomes weaker, accordingly it needs smaller temperature to buckle.

Figure 3.

Effects of the length-to-depth ratio $\gamma =L/\tilde{H}$ and external electric voltage $E_0=V_a/({\tilde{H}}^2{E}_{\mathrm{pz}})$ on the critical buckling temperature $T_{\mathrm{cr}}$ of FGGPL strengthened piezoelectric beams for four different patterns: (a) pattern 1, (b) pattern 2, (c) pattern 3, and (d) pattern 4 ( $F^{\mathrm{GP}}=0.1,e_0=300\times {10}^3,\mathrm{\Delta }C=0.01$ ).

Figure 4 represents the influences of the length-to-depth ratio $\gamma$ and moisture concentration $\mathrm{\Delta }C$ on the critical buckling temperature $T_{\mathrm{cr}}$ of FGGPL strengthened piezoelectric beams where the four patterns (a) pattern 1, (b) pattern 2, (c) pattern 3, and (d) pattern 4 are considered. It is observed that the critical buckling temperature decreases as both of the ratio $\gamma$ and moisture concentration $\mathrm{\Delta }C$ increase. Furthermore, the effects of $\mathrm{\Delta }C$ on $T_{\mathrm{cr}}$ become more pronounced for larger values of the ratio $\gamma$ . Indeed, increasing the humidity weakens the beam stiffness which leads to less buckling temperature.

Figure 4.

Effects of the length-to-depth ratio $\gamma$ and moisture concentration $\mathrm{\Delta }C$ on the critical buckling temperature $T_{\mathrm{cr}}$ of FGGPL strengthened piezoelectric beams for four different patterns (a) pattern 1, (b) pattern 2, (c) pattern 3, and (d) pattern 4 ( $F^{\mathrm{GP}}=0.1,e_0=300\times {10}^3,E_0=0.01$ ).

Figure 5 illustrates the changes of the critical buckling temperature $T_{\mathrm{cr}}$ of FGGPL strengthened piezoelectric beams regarding different values of graphene volume fraction $F^{\mathrm{GP}}$ against the length-to-depth ratio $\gamma$ for four different patterns 1, 2, 3, and 4 in the corresponding figures a, b, c, and d, respectively. It is noticed that the beam stiffness increases as the graphene volume fraction $F^{\mathrm{GP}}$ increases leading to an increment in the buckling temperature $T_{\mathrm{cr}}$ . It is obvious for pattern 1 that the impacts of $F^{\mathrm{GP}}$ on $T_{\mathrm{cr}}$ become more noticeable for smaller length-to-depth ratios $\gamma$ .

Figure 5.

Effects of the length-to-depth ratio $\gamma$ and graphene volume fraction $F^{\mathrm{GP}}$ on the critical buckling temperature $T_{\mathrm{cr}}$ of FGGPL strengthened piezoelectric beams for four different patterns: (a) pattern 1, (b) pattern 2, (c) pattern 3, and (d) pattern 4 ( $\mathrm{\Delta }C=0.01,e_0=300\times {10}^3,E_0=0.01$ ).

Figure 6 displays the effects of the external electric voltage $E_0$ and piezoelectric multiple $e_0$ on the critical buckling temperature $T_{\mathrm{cr}}$ of FGGPL strengthened piezoelectric beams for four different patterns (a) pattern 1, (b) pattern 2, (c) pattern 3, and (d) pattern 4. It is obvious for all patterns that the critical buckling temperatures are linearly varied and increase as both $E_0$ and $e_0$ increase. Moreover, the impacts of piezoelectric multiple $e_0$ on $T_{\mathrm{cr}}$ are more evident for higher values of $E_0$ .

Figure 6.

Effects of the external electric voltage $E_0$ and piezoelectric multiple $e_0$ on the critical buckling temperature $T_{\mathrm{cr}}$ of FGGPL strengthened piezoelectric beams for four different patterns: (a) pattern 1, (b) pattern 2, (c) pattern 3, and (d) pattern 4 ( $\mathrm{\Delta }C=0.01,\gamma =10,F^{\mathrm{GP}}=0.01$ ).

Conclusions

This research is devoted to investigate the electro-thermal buckling behavior of piezoelectric FG graphene platelets beam (PFGGPL) which subjected to an external electric voltage and moisture conditions. In accordance with Halpin-Tsai model, the effective Young’s modulus of the beam was predestined. Depending on the rule of the mixture, Poisson’s ratio, mass density, and piezoelectric properties were measured under four FGGPL distribution types. The displacement components were expressed by using a refined two-unknown beam theory taking into account the shear deformation as well as thickness stretching effect. Moreover, the stability differential equations were derived by employing the principle of virtual work containing both of the thermal, moisture, and electric loads. Furthermore, the present formulations were validated by comparing the related outcomes with those of other published articles. It was concluded for all FGGPL distribution types that the critical buckling temperature increases by increasing each of the external electric voltage, the piezoelectric multiple, and the graphene volume fraction. However, the critical buckling temperature decreases by increasing both of the length-to-depth ratio and the moisture concentration.