Influence of Applied Loads on Free Vibrations of Functionally Graded Material Plate-Shell Panels

1. Introduction

Plate/shells panels have applications in a wide range of engineering situations, such as aeronautics and aerospace structural systems, among many others. Structures made of composite materials have been widely used to satisfy high performance demands. In such structures, stress singularities may occur at the interface between two different materials. In functionally graded structures the smooth and continuous variation of the properties, from one surface to the other, eliminates abrupt changes in the stress and displacement distributions. In addition, ceramic phase with the low thermal conductivity can resist high thermal environment while metal phase is strong with mechanical load.

Thousands of research works in plate/shells panels composite structures had been published in the last five decades. More recently also a great number of works had been published in functionally graded material structures with the same type of geometry. For this reason, only some representative works related with the objective of the present investigation are mentioned dealing with functionally graded panels structures [1,2,3,4,5,6,7,8,9,10,11,12,13].

In this paper it is presented the analysis of functionally graded material plate/shells panels subjected to general loadings. The formulation is developed for the calculation of natural frequencies and to study the influence of the applied loads on these natural frequencies. Thus, the formulation has to allow for different analyses: geometrically nonlinear deformation with the through-thickness stress distribution calculations, linear buckling and free vibrations. A finite element model based on a higher order shear deformation theory and using a non-conforming triangular flat plate/shell element with three nodes, and eight degrees of freedom per node is adopted to study the plate/shell structures. This finite element model reveals to be very efficient in the analysis of the proposed structures – easy to model the geometry, low number of total degrees of freedom, good accuracy of results.

An in-house computer program was developed, and the solutions of some illustrative examples are performed and the results are presented and discussed.

2. Formulation of the P-FGM Model

An FGM is made by mixing two distinct isotropic material phases, for example a ceramic and a metal. In this work the material properties of an FGM plate/shell structure are assumed to change continuously throughout the thickness (Figure 1), according to the volume fraction of the constituent materials, given by the Power-Law function – Bao and Wang [14]. In addition, the continuous variation of the materials mixture is approximated by using a certain number of virtual layers k throughout the thickness direction - layer approach. The volume fraction of the ceramic and metal phases for each virtual layer is defined according to the power-law:

$$V_c^k={\left(0.5+\frac{\bar{z}}{h}\right)}^{ p}; V_m^k=1.0-V_c^k$$

(1)

where $\bar{z}$ is the thickness coordinate of mid-surface of each layer.

Once the volume fractions of the ceramic and the metal, $V_c^k$ and $V_m^k$, are defined, the Young’s modulus E, the shear modulus $\mathrm{G}$, the Poisson ratio $\nu$, or the mass density $\rho$ of each virtual layer k of the FGM plate/shell can be determined:

$$\begin{gathered} E_k=V_c^k{\mathrm{E}}_c+V_m^k{\mathrm{E}}_m ; {\mathrm{G}}_k=V_c^k{\mathrm{G}}_c+V_m^k{\mathrm{G}}_m \\ {\nu }_k=V_c^k {\nu }_c+V_m^k {\nu }_m ; {\rho }_k=V_c^k {\rho }_c+V_m^k {\rho }_m \end{gathered}$$

(2)

2.1. Displacement Fields and Strain-Displacement Relations

The displacement field of the plate/shell domain is based on the Reddy’s third-order shear deformation theory, [12,13]:

$$\begin{gathered} u\left(x,y,z\right)=u_0\left(x,y\right)-z {\theta }_y\left(x,y\right)+z^3{\mathrm{c}}_1\left[{\theta }_y\left(x,y\right)-\frac{\partial {\mathrm{w}}_0}{\partial \mathrm{x}}\right] \\ v\left(x,y,z\right)=v_0\left(x,y\right)+z {\theta }_x\left(x,y\right)+z^3{\mathrm{c}}_1\left[-{\theta }_x\left(x,y\right)-\frac{\partial {\mathrm{w}}_0}{\partial \mathrm{y}}\right] \\ w\left(x,y,z\right)=w_0\left(x,y\right) \end{gathered}$$

(3)

where $u_0, v_0, w_0$ are displacements of a generic point in the middle plane referred to the local axes - x,y,z directions, ${\theta }_x, {\theta }_y$ are the rotations of the normal to the middle plane, about the $x$ axis (clockwise) and $y$ axis (anticlockwise), $\frac{\partial w_0}{\partial x} , \frac{\partial w_0}{\partial y}$ are the slopes of the tangents of the deformed mid-surface in x,y directions, and $c_1={4/3 h}^2$, with $h$ denoting the total thickness of the plate/shell structure.

For the geometrically nonlinear behaviour analysis, the strain tensor is conveniently represented in terms of the linear and nonlinear parts of the strain tensor as:

$$\epsilon =\left\{ {\epsilon }^L+{\epsilon }^{NL}\right\}$$

(4)

The linear strain components associated with the displacement fields defined above can be written in the form:

$${\epsilon }^L=\left\{\begin{array}{c}{\mathsf{\epsilon }}^m+z{\mathsf{\epsilon }}^b+z^3 {\mathsf{\epsilon }}^{\ast b} \\ {\mathsf{\epsilon }}^s+z^{2 }{\mathsf{\epsilon }}^{\ast s}\end{array}\right\}$$

(5)

Therefore, we can write in the detailed form:

$$\begin{gathered} {\epsilon }_{xx}=\frac{\partial u_0}{ \partial x}-z \frac{ \partial {\theta }_y}{ \partial \mathrm{x}}+z^3{c}_1\left(\frac{\partial {\theta }_y}{\partial x}-\frac{{\partial }^2{w}_0}{\partial x^2}\right)+\frac{1}{2}\left[{\left(\frac{\partial u_0}{\partial x}\right)}^2+{\left(\frac{\partial v_0}{\partial x}\right)}^2+{\left(\frac{\partial w_0}{\partial x}\right)}^2\right] \\ {\epsilon }_{yy}=\frac{\partial v_0}{\partial y}+z\frac{ \partial {\theta }_x}{\partial y}+z^3{c}_1\left(-\frac{\partial {\theta }_x}{\partial y}-\frac{{\partial }^2{w}_0}{\partial y^2}\right)+\frac{1}{2}\left[{\left(\frac{\partial u_0}{\partial y}\right)}^2+{\left(\frac{\partial v_0}{\partial y}\right)}^2+{\left(\frac{\partial w_0}{\partial y}\right)}^2\right] \\ {\gamma }_{xy}=\left(\frac{\partial u_0}{\partial y}+\frac{\partial v_0}{\partial x}\right)+z \left(-\frac{\partial {\theta }_y}{\partial y}+\frac{\partial {\theta }_x}{\partial x}\right)+z^3{c}_1\left( \frac{\partial {\theta }_y}{\partial y}-\frac{\partial {\theta }_x}{\partial x}-2\frac{{\partial }^2{w}_0}{\partial x \partial y}\right)+ \\ \left[\left(\frac{\partial u_0}{\partial x}\frac{\partial u_0}{\partial y}\right)+\left(\frac{\partial v_0}{\partial x}\frac{\partial v_0}{\partial y}\right)+\left(\frac{\partial w_0}{\partial x}\frac{\partial w_0}{\partial y}\right)\right]~ \\ {\gamma }_{xz}=-{\theta }_y+\frac{\partial w_0}{\partial x}+z^2{c}_2\left({\theta }_y-\frac{\partial w_0}{\partial x}\right) \\ {\gamma }_{yz}={\theta }_x+\frac{\partial w_0}{\partial y}+z^2{c}_2\left(-{\theta }_x-\frac{\partial w_0}{\partial y}\right) \end{gathered}$$

(6)

2.2. Constitutive Relations for the FGM Structures

The stress-strain relations are defined based on the stress-strain relations in each layer k, which can be written as follows:

$${\sigma }_k=Q_k {\epsilon }_k$$

(7)

where ${\sigma }_k={\left\{\begin{array}{ccccc}{\sigma }_x& {\sigma }_y& {\sigma }_{xy}& {\tau }_{xz}& {\tau }_{yz}\end{array}\right\}}^T$ is the stress vector, ${\epsilon }_k={\left\{\begin{array}{ccccc}{\epsilon }_x& {\epsilon }_y& {\gamma }_{xy}& {\gamma }_{xz}& {\gamma }_{yz}\end{array}\right\}}^T$ is the strain vector and $Q_k$ is the elasticity matrix of the k layer:

$$Q_k={{\left[\begin{array}{ccccc}{Q}_{11}& Q_{12}& 0& 0& 0\\ Q_{21}& Q_{22}& 0& 0& 0\\ 0& 0& Q_{66}& 0& 0\\ 0& 0& 0& Q_{44}& 0\\ 0& 0& 0& 0& Q_{55}\end{array}\right]}_k}_{ }$$

(8)

$$Q_{11_k}={Q_{22}}_k=\frac{E_k}{1-{\nu }_k^2};Q_{12_k}={Q_{21}}_k=\frac{{\nu }_k E_k}{1-{\nu }_k^2};$$

$$Q_{44_k}=Q_{55_k}=Q_{66_k}=\frac{E_k}{2 (1+\nu )}=G_k$$

(9)

By integrating through the thickness, the linear elastic constitutive equation for each virtual layer k, can be written as:

$${\widehat{\sigma }}_k={\widehat{D}}_k {\epsilon }_k={\left\{\begin{array}{c}N\\ M\\ M^{\ast }\\ Q\\ Q^{\ast }\end{array}\right\}}_k={\left[\begin{array}{ccccc}A& B& D& 0& 0\\ B& C& E& 0& 0\\ D& E& G& 0& 0\\ 0& 0& 0& A_s& C_s\\ 0& 0& 0& C_s& E_s\end{array}\right]}_k{\left\{\begin{array}{c}{\epsilon }^m\\ {\epsilon }^b\\ {\epsilon }^{\ast b}\\ {\epsilon }^s\\ {\epsilon }^{\ast s}\end{array}\right\}}_k$$

(10)

where ${\widehat{\sigma }}_k$are the resultant forces and moments, ${\widehat{D}}_k \mathrm{i}s$the constitutive matrix, and the sub-matrices of $\widehat{D}$ are given by:

$$\begin{gathered} {\left(A,B,C,D,E,G\right)}_k={\bar{Q}}_k {\int }_{z_{k-1}}^{z_k}(1,z,z^2,z^3,z^4,z^6) dz; \\ {\left(A_s,C_s,E_s\right)}_k={\bar{Q}}_k {\int }_{z_{k-1}}^{z_k}(1,z^2,z^4) dz. \end{gathered}$$

(11)

3. Finite Element Approach

A non-conforming triangular plate/shell finite element model having three nodes and eight degrees of freedom per node is used. The degrees of freedom at each node are the displacements, the slopes and the rotations, $d_{\mathrm{i}}={\left\{u_0 v_0 w_0 -\frac{\partial w_0}{\partial y} \frac{\partial w_0}{\partial x} {\theta }_x {\theta }_y {\theta }_z\right\}}_{\mathrm{i}}^T.$

The rotation ${\theta }_{z\mathrm{i}}$ is introduced to consider a fictitious stiffness coefficient $K_{\theta z}$ to eliminate the problem of a singular stiffness matrix for general shape structures [15]. The finite element formulation is obtained by the superposition of a membrane element with a bending element, and the side 1-2 of the triangle is coincident with the axis x of local system attached to the finite element, with origin at node 1 (Figure 2).

The element displacement vector at any point of mid-surface in local coordinates (x,y) are obtained in terms of the corresponding nodal variables through linear shape functions $L_{\mathrm{i}}$, and the transverse displacement and slopes are expressed in terms of corresponding nodal variables through cubic shape functions ${{}_{j}N}_{\mathrm{i}}$, [15], which are also given in terms of area co-ordinates [15].

The displacement field can be represented in matrix form as:

$$u=Z \left({\sum }_{\mathrm{i}=1}^3{N}_{\mathrm{i}} d_{\mathrm{i}}\right)=\mathrm{Z} \mathrm{N} {\mathrm{a}}^e ; a^e={\left\{\begin{array}{ccc}{d}_1& d_2& d_3\end{array}\right\}}^T$$

(12)

where the matrix of shape functions and the appropriate matrix Z containing powers of $z_k$, are:

$$N_{\mathrm{i}}=\left[\begin{array}{cccccccc}{L}_{\mathrm{i}}& 0& 0& 0& 0& 0& 0& 0\\ 0& L_{\mathrm{i}}& 0& 0& 0& 0& 0& 0\\ 0& 0& {{}_{1}N}_{\mathrm{i}}& {{}_{2}N}_{\mathrm{i}}& {{}_{3}N}_{\mathrm{i}}& 0& 0& 0\\ 0& 0& -\frac{\partial {{}_{1}N}_{\mathrm{i}}}{\partial y}& -\frac{{\partial }_2{N}_{\mathrm{i}}}{\partial y}& -\frac{\partial {{}_{3}N}_{\mathrm{i}}}{\partial y}& 0& 0& 0\\ 0& 0& \frac{\partial {{}_{1}N}_{\mathrm{i}}}{\partial x}& \frac{\partial {{}_{2}N}_{\mathrm{i}}}{\partial x}& \frac{\partial {{}_{3}N}_{\mathrm{i}}}{\partial x}& 0& 0& 0\\ 0& 0& 0& 0& 0& L_{\mathrm{i}}& 0& 0\\ 0& 0& 0& 0& 0& 0& L_{\mathrm{i}}& 0\end{array}\right]$$

(13)

$$Z=\left[\begin{array}{cccccccc}1& 0& 0& 0& -z^3{c}_1& 0& -z+z^3{c}_1& 0\\ 0& 1& 0& z^3{c}_1& 0& z-z^3{c}_1& 0& 0\\ 0& 0& 1& 0& 0& 0& 0& 0\end{array}\right]$$

(14)

4. Virtual Work Principle

4.1. Static Analysis

The governing equations of the nonlinear problem are obtained from the virtual work principle in conjugation with an updated Lagrangian formulation [16,17]. A reference configuration is associated with time t and the updated configuration is associated with the current time ${\mathrm{t}}^{\prime }=\mathrm{t}+\Delta \mathrm{t}$.

$$\begin{gathered} {\sum }_{\mathrm{k}=1}^{\mathrm{N}}\left\{{\int }_{{\mathrm{t}}_{{\mathrm{A}}^{\mathrm{e}}}}{\int }_{{\mathrm{h}}_{\mathrm{k}-1}}^{{\mathrm{h}}_{\mathrm{k}}}\mathsf{\delta }\left({}_{\mathrm{t}}{\mathsf{\epsilon }}_{\mathrm{k}}^{\mathrm{L}}\right){\mathrm{Q}}_{\mathrm{k}}{}_{\mathrm{t}}\mathsf{\epsilon }_{\mathrm{k}}^{\mathrm{L}}\mathrm{d}\mathrm{z}{}^{\mathrm{t}}\mathrm{d}\mathrm{A}^{\mathrm{e}}+{\int }_{{\mathrm{t}}_{{\mathrm{A}}^{\mathrm{e}}}}{\int }_{{\mathrm{h}}_{\mathrm{k}-1}}^{{\mathrm{h}}_{\mathrm{k}}}\mathsf{\delta }\left({}_{\mathrm{t}}{\mathsf{\epsilon }}_{\mathrm{k}}^{\mathrm{N}\mathrm{L}}\right){}^{\mathrm{t}}\mathsf{\sigma }_{\mathrm{k}}^{\mathrm{L}}\mathrm{d}\mathrm{z}{}^{\mathrm{t}}\mathrm{d}\mathrm{A}^{\mathrm{e}}\right\}={}_{\mathrm{t}}{}^{\mathrm{t}+∆\mathrm{t}}{\mathsf{\Pi }}_{\mathrm{e}}- \\ {\sum }_{\mathrm{k}=1}^{\mathrm{N}}{\int }_{{\mathrm{t}}_{{\mathrm{A}}^{\mathrm{e}}}}{\int }_{{\mathrm{h}}_{\mathrm{k}-1}}^{{\mathrm{h}}_{\mathrm{k}}}\mathsf{\delta }\left({}_{\mathrm{t}}{\mathsf{\epsilon }}_{\mathrm{k}}^{\mathrm{L}}\right){}^{\mathrm{t}}\mathsf{\sigma }_{\mathrm{k}}^{\mathrm{L}}\mathrm{d}\mathrm{z}{}^{\mathrm{t}}\mathrm{d}\mathrm{A}^{\mathrm{e}} \end{gathered}$$

(15)

where ${}_{\mathrm{t}}{}^{\mathrm{t}+∆\mathrm{t}}{\Pi }_e$ is the element external virtual work.

The membrane, bending and shear strains, as well the higher order bending, and shear strains can be represented by:

$${\mathsf{\epsilon }}_{\mathrm{m}}={\mathbf{B}}^m{\mathbf{a}}^e ; {\mathsf{\epsilon }}_{\mathrm{b}}={\mathbf{B}}^b\mathbf{z}{\mathbf{a}}^e ; {\mathsf{\epsilon }}_{\mathrm{s}}={\mathbf{B}}^s{\mathrm{a}}^e; {\mathsf{\epsilon }}_{\mathrm{b}}^{\mathrm{*}}={\mathbf{B}}^{\mathbf{*}b}{\mathrm{z}}^3{\mathrm{a}}^e ; {\mathsf{\epsilon }}_{\mathrm{s}}^{\mathrm{*}}={\mathrm{B}}^{\mathrm{*}s}{\mathrm{z}}^2{\mathbf{a}}^e$$

(16)

where the${ \mathbf{B}}^m$, ${\mathbf{B}}^b$, ${\mathbf{B}}^{\ast b}$, ${\mathbf{B}}^s$, ${\mathbf{B}}^{\ast s}$, are components of the strain-displacement matrix $\mathbf{B},$ and are given explicitly in Moita et al. [18].

The virtual work principle applied to a finite element, can be written in the following matrix form:

$${\int }_{{\mathrm{t}}_{{\mathrm{A}}^{\mathrm{e}}}}{\mathbf{B}}^{\mathit{T}}\widehat{\mathbf{D}}\mathbf{B}{∆\mathrm{a}}^e{{}^{\mathrm{t}}\mathrm{d}\mathrm{A}}^{\mathrm{e}}+{\int }_{{\mathrm{t}}_{{\mathrm{A}}^{\mathrm{e}}}}{\mathbf{G}}^{\mathit{T}}\widehat{\mathsf{\tau }}\mathbf{G} ∆{\mathrm{a}}^e{{}^{\mathrm{t}}\mathrm{d}\mathrm{A}}^{\mathrm{e}}={\int }_{A_e}{\mathbf{N}}_T\mathbf{p}{dA}_e-{\int }_{{\mathrm{t}}_{{\mathrm{A}}^{\mathrm{e}}}} {\mathbf{B}}^T{}^{\mathrm{t}}\widehat{\mathsf{\sigma }} {{}_{\mathrm{t}}\mathrm{d}\mathrm{A}}^{\mathrm{e}}$$

(17)

From this equation, it can be obtained the element linear stiffness matrix ${\mathbf{K}}_L^e,$ the element geometric stiffness matrix ${\mathbf{K}}_{\sigma }^e$, the element external load vector ${\mathbf{F}}_{ex\mathrm{t}}^e$ and the element internal force vector ${\mathbf{F}}_{\mathrm{i}n\mathrm{t}}^e$ [18].

The element geometric stiffness matrix${\mathbf{K}}_{\sigma }^e$ is given by:

$${\mathbf{K}}_{\mathsf{\sigma }}^{\mathrm{e}}={\int }_{{\mathrm{t}}_{{\mathrm{A}}^{\mathrm{e}}}}\left({\mathbf{G}}^{\mathrm{T}}\widehat{\mathsf{\tau }}\mathbf{G}\right){{}^{\mathrm{t}}\mathrm{dA}}^{\mathrm{e}}$$

(18)

$${\mathrm{G}}_{\mathrm{i}}=\left[\begin{array}{cccccccc}\frac{\partial L_{\mathrm{i}}}{\partial x}& 0& 0& 0& 0& 0& 0& 0\\ \frac{\partial L_{\mathrm{i}}}{\partial y}& 0& 0& 0& 0& 0& 0& 0\\ 0& \frac{\partial L_{\mathrm{i}}}{\partial x}& 0& 0& 0& 0& 0& 0\\ 0& \frac{\partial L_{\mathrm{i}}}{\partial y}& 0& 0& 0& 0& 0& 0\\ 0& 0& \frac{\partial {}_{{}_{1}N_{\mathrm{i}}}}{\partial x}& \frac{\partial {}_{{}_{2}N_{\mathrm{i}}}}{\partial x}& \frac{\partial {}_{{}_{3}N_{\mathrm{i}}}}{\partial x}& 0& 0& 0\\ 0& 0& \frac{\partial {}_{{}_{1}N_{\mathrm{i}}}}{\partial y}& \frac{\partial {}_{{}_{2}N_{\mathrm{i}}}}{\partial y}& \frac{\partial {}_{{}_{3}N_{\mathrm{i}}}}{\partial y}& 0& 0& 0\end{array}\right]$$

(19)

$$\widehat{\tau }=\left[\begin{array}{cccccc}{N}_{xx}& N_{xy}& 0& 0& 0& 0\\ N_{xy}& N_{yy}& 0& 0& 0& 0\\ 0& 0& N_{xx}& N_{xy}& 0& 0\\ 0& 0& N_{xy}& N_{yy}& 0& 0\\ 0& 0& 0& 0& N_{xx}& N_{xy}\\ 0& 0& 0& 0& N_{xy}& N_{yy}\end{array}\right]$$

(20)

To solve for general structures, these matrices and vectors, initially computed in the local coordinate system attached to the element, are transformed to global coordinate system in the usual way [15]. Afterwards, by adding the contributions of all the elements in the domain, the system equilibrium equations can be obtained as:

$${}_{\mathrm{t}+\Delta \mathrm{t}}{}^{\mathrm{t}+\Delta \mathrm{t}}\left({\mathbf{K}}_L+{\mathbf{K}}_{\sigma }\right){ }^{\mathrm{i}-1} \left(\Delta \mathbf{q}\right){ }^{\mathrm{i}}={}^{\mathrm{t}+\Delta \mathrm{t}}\mathbf{F}_{ex\mathrm{t}}-{}_{\mathrm{t}+\Delta \mathrm{t}}{}^{\mathrm{t}+\Delta \mathrm{t}}\left({\mathbf{F}}_{\mathrm{int}}\right){ }^{\mathrm{i}-1}$$

(21)

Using the Newton-Raphson incremental-iterative method [20], the incremental equilibrium path is obtained, and in the case of snap-through occurrence the automatic arc-length method is used [19] being the Eq. (21) written in the form:

$${}_{\mathrm{t}+\Delta \mathrm{t}}{}^{\mathrm{t}+\Delta \mathrm{t}}\left({\mathbf{K}}_L+{\mathbf{K}}_{\sigma }\right){ }^{\mathrm{i}-1} \left(\Delta \mathbf{q}\right){ }^{\mathrm{i}}=\left({\lambda }^{\mathrm{i}-1}+\Delta {\lambda }^{\mathrm{i}}\right){\mathbf{F}}_{ex\mathrm{t}}^0-{}_{\mathrm{t}+\Delta \mathrm{t}}{}^{\mathrm{t}+\Delta \mathrm{t}}\left({\mathbf{F}}_{\mathrm{int}}\right){ }^{\mathrm{i}-1}$$

(22)

and an additional equation is employed to constrain the length of a load step:

$${{\left(\mathsf{\Delta }\mathbf{q}\right)}^{\mathrm{i}}}^{\mathrm{T}}{\left(\mathsf{\Delta }\mathbf{q}\right)}^{\mathrm{i}}={\left(\mathsf{\Delta }\mathcal{l}\right)}^2$$

(23)

where ${\mathbf{F}}_{ex\mathrm{t}}^0$ is a fixed (reference) load vector, $\mathsf{\lambda }$ is a load factor, $\mathsf{\Delta }\mathsf{\lambda }$ is the incremental load factor within the load step, and $\mathsf{\Delta }\mathcal{l}$ is the arc-length.

4.2. Linear Buckling Analysis

For linear mechanical buckling analysis, only one load increment is used, and equation (21) becomes:

$$\left({\mathbf{K}}_{\mathrm{L}}^0+{\mathbf{K}}_{\mathsf{\sigma }}^0\right)\mathbf{q}={\mathbf{F}}_{\mathrm{e}\mathrm{x}\mathrm{t}}^0$$

(24)

The linear elastic buckling standard eigenvalue problem is then given by:

$$\left({\mathbf{K}}_L\hspace{0.33em}+\hspace{0.33em}\hspace{0.33em}{\mathsf{\lambda }}_{\mathrm{c}\mathrm{r}}\hspace{0.33em}{\mathbf{K}}_{\mathsf{\sigma }}^0\right)\mathbf{q}\hspace{0.33em}=0 , {}_{cr}{\mathbf{F}}_{ex\mathrm{t}}={\lambda }_{cr} {\mathbf{F}}_{ex\mathrm{t}}^0$$

(25)

4.3. Free Vibration Analysis

The Hamilton principle is used to formulate the governing equations for vibration analysis:

$${\int }_{t_1}^{t_2}\left[{\sum }_{k=1}^N\left({\int }_{A^e} {\int }_{h_{k-1}}^{h_k} \mathsf{\delta }\left({\bar{\epsilon }}_k^L\right) {\bar{Q}}_k {\bar{\epsilon }}_k^L dz {}_{t}d{A}^e-\right.\right.\left.\left.{\int }_{A^e} {\int }_{h_{k-1}}^{h_k} \delta {\dot{u}}^T {\rho }_k \dot{u} dz {}_{t}d{A}^e \right)\right]dt=0$$

(26)

$${\int }_{t_1}^{t_2}\left[{\sum }_{k=1}^N \left({\int }_{A^e}{\int }_{h_{k-1}}^{h_k}\delta {{\mathrm{a}}^e}^T{B}^T\bar{Q} B a^e dz d{A}^e- {\int }_A\delta {{\dot{\mathrm{a}}}^e}^T{N}^{T }{\rho }_k\left({\int }_{h_{k-1}}^{h_k}{Z}^T\mathrm{Z} \mathrm{d}\mathrm{z}\right) N{\dot{\mathrm{a}}}^e d{A}^e\right)\right] dt=0$$

(27)

where the element mass matrix M^e is given by:

$$M^e={\int }_A{N}^{\mathrm{T} }{\sum }_{k=1}^N{\rho }_k\left({\int }_{h_{k-1}}^{h_k}{Z}^T\mathrm{Z} \mathrm{d}\mathrm{z}\right) N\mathrm{d}{\mathrm{A}}^e$$

(28)

After assembling all the elements, the natural frequencies and the respective mode shapes are obtained by solving the standard eigenvalue problem:

$$\left(\left[{\mathbf{K}}_L\right]-{\omega }^2 \left[\mathbf{M}\right]\right)\mathbf{q}=0$$

(29)

The free vibrations of loaded panels are obtained considering the corresponding stiffness matrix by the following equation:

$$\left({}_{\mathrm{t}}{}^{\mathrm{t}}\left[{\mathbf{K}}_L+{\mathbf{K}}_{\sigma }\right]-{\omega }^2{ }_{\mathrm{t}}^{\mathrm{t}}\left[\mathbf{M}\right]\right)\mathbf{q}=0$$

(30)

5. Applications

A set of numerical aplications concerning the objective of the present work are presented. With this purpose we present also a few numerical applications, considering different type of analyses whose are necessary to the calculations of the final objective.

5.1. Nonlinear Analysis of a Simply Supported and Clamped FGM Square Plate under Uniform Pressure

Simply supported (SS) and clamped (CL) FGM square plates subjected to a uniform pressure load are considered in this first numerical application. The geometry is defined by the side dimension a=1.0 m and the thickness h = 0.01 m. Three gradient indexes are considered: p=0.5, p=1 and p=2. The plate, discretized through the thikcness by using 20 virtual layers, is made of constituents zirconia and aluminium (E_c = 151x10⁹ Pa, ${\nu }_c$= 0.3, ρ_c=5700 kg/m³, E_m =70 x 10⁹ Pa, ${\nu }_m$=0.3, ρ_m=2702 kg/m³). An incremental-iterative loading is applied and the load-displacement paths, obtained for the centre of the plate, are presented in Figure 3. This analysis had been validated in Moita et al. [13], although for different geometry and material properties. It is observed that the displacements increase when the gradient index increases. It happens because E_c>E_m, and as p increases the volume fraction of ceramic decreases.

For the last load increment, considering the simply supported plate with gradient index p=1, the through-the-thickness distribution of total normal and transverse shear stresses at Gauss point nearest of centre point of the plate are shown in Figure 4 and Figure 5. It is observed a stress distribution of stresses without abrupt changes in the through-the-thickness as expected.

5.2. Nonlinear Analysis of a Hinged and Clamped FGM Cylindrical Panels under a Centre point Load

Two FGM cylindrical shell panels (Figure 6), one with the straight sides simply supported (hinged) and the curved sides free, and the other with all sides clamped, are analyzed in this subsection. The geometry is defined by: R= 2.54 m, L = 0.508 m and subtended angle 2$\theta$=0.2 rad, H=0.0127 m, and it is modelled by 392 triangular finite elements. The constituents are ZrO₂ and Ti–6Al–4V with the mechanical properties, respectively: E_c = 168.0 x 109 N/m²,${\nu }_c$= 0.298, ρ_c=5700 kg/m³ , E_m = 105.7 x109 N/m²,${\nu }_m$= 0.298, ρ_m=4429 kg/m³. For this case of geometry and boundary conditions of the cylindrical panel, the snap-through can occur. The results obtained using the present model for the transverse displacement paths, are shown in Figure 7a,b and Figure 8a,b, respectively for the hinged cylindrical panel and for the clamped cylindrical panel. These analysis have been validated in Moita et al. [13], although for different geometry and material properties. It is observed from the paths obtained for the hinged cylindrical panel, that the snap-through do occur, while for the clamped FGM panels, the snap-through do not occur. Also, as expected, the displacements increase when the gradient index increase, considering that Ec>Em. Considering the same cylindrical panel with both boundary conditions, but under radial pressure load, the corresponding load-displacement paths are shown in Figure 7b and Figure 8b.

5.3. Critical Loads of a Clamped FGM Square Plate

The linear buckling loads are obtained first in order to know which values should be considered for the calculations of the free vibrations under uniaxial in-plane applied load.

For the purpose of comparison, a square (a×a) FGM plate with all edges movable and clamped, made of a mixture of zirconia and aluminium (Al/ZrO2), is considered. The material properties for the zirconia and aluminium are respectively, Ec=151.0×109 N/m2,${\nu }_c=0.3$, Em = 70.0×109 N/m2, ${\nu }_m=0.3$, side dimension are a=0.2 m and thickness H=0.01m. The FGM plates are under uniaxial compression and the power law index-p ranges from full ceramic to metal. The critical buckling parameters, defined as ${\lambda }_{cr}=N_{xx}{a}^2/D_0$, with $D_0=E_m{h}^3/12(1-{\nu }^2)$, are shown in Table 1. From this table it is observed an excellent agreement between the results obtained with the present model (PM) compared with those obtained by Nguyen-Xuan et al. [20] using the MITC4 model, for the first mode. For benchmarking purposes, the eigenvalues of some higher modes are also presented. Again, as Ec>Em, there is a decreasing in the buckling loads when the p-index increase.

Also for comparison purposes, it is considered a square FGM plate with all edges movable and clamped, but now with H=0.01 m and sides with length a=1.0 m. The material propoetirs are the following: Ec=380.0×109 N/m2,${\nu }_c=0.3$, Em = 70.0×109 N/m2, ${\nu }_m=0.3$. The FGM plates are under uniaxial compression and the power law index-p ranges from full ceramic to metal. The critical buckling parameter is now defined as ${\lambda }_{cr}=N_{xx}{a}^2/E_c{H}^3$, and the results obtained with the present model (PM) are compared in Table 2 with the alternative solutions, for the first mode, obtained analytically by Wu et al. [21] using the first-order shear deformation theory.

5.4. Free Vibrations of a Simply Supported FGM Square Plate

The free vibrations of a simply supported FGM square plate are now investigated. The material properties for the zirconia and aluminium are respectively, the following: E_c=151.0×109 N/m²,${\nu }_c=0.3$, ρ_c=5700 kg/m³, and E_m = 70.0×109 N/m², ${\nu }_m=0.3$, ρ_m=2702 kg/m³. The side dimension is a=0.05 m and the thickness H=0.01m. The non-dimensional frequencies $\bar{\omega }=\omega h\sqrt{{\rho }_m/E_m}$obtained with the present model (PM) are presented in Table 3 and compared with the alternative solutions obtained by Loc et al [9] using a NURBS-based isogeometric approach and a HSDT model. For pure ceramic and pure metal the results are very similar. The discrepancies between the results for FGM plates are attributed to the different homogenization considered - Loc et al. [9] follow the Mori–Tanaka scheme, while the present model uses the rule of mixtures.

The next investigation is to achieve how the free vibrations of a FGM plate with all edges clamped and power-law index p=1.0, are affected by two types of loading, the in-plane loading and transverse loading. The material properties are: E_c=380.0×109 N/m²,${\nu }_c=0.3$, ρ_m=3000 kg/m³, E_m = 70.0×109 N/m², ${\nu }_m=0.3$, ρ_m=2702 kg/m³. The geometry is defined by a=1.0 m and H=0.01 m. The results obtained are given in Table 4a–c. First, the free vibrations are calculated for the unloaded plate, and afterwards for the three cases of loading:

a) uniaxial in-plane uniform load from p_x=0 to a final p_x=2.0; 4.0; 10.0 MPa, and the geometrically nonlinear behavior is considered;

b) transverse uniform load from p_z=0 to a final p_z=20.0 ; 40.0; 60.0 kPa;

c) uniaxial in-plane uniform load combined with transverse uniform load, ranging from zero to the previously defined magnitudes.

From the Table 4a–c, it is observed that the uniform compressive loading decreases the natural frequencies while the uniform transverse pressure increases the natural frequencies. In addition, it is observed when both types of loads are applied the first frequency increases, but less than when only the transverse pressure is applied, and the other frequencies increase or decrease depending of the values of two types of loadings. Also, it is shown that the magnitude of decreasing/increasing is directly connected with the increase of loading.

Considering the gradient index p=1.0, the influence on the fundamental frequency of the applied loadings is shown graphically in Figure 9a,b for the case of clamped plate, according to Table 4a–c.

5.5. Critical Loads of Hinged and Clamped FGM Cylindrical Panels

To validate the present model, the buckling analysis of a square FGM panel made of a mixture of aluminum and zirconia with the following properties E_m = 70.0 x109 N/m², ${\nu }_m$= 0.3, E_c = 151.0 x 109 N/m², ${\nu }_c$= 0.3, is considered. The geometry of the panel is length L = 0.1m, angle 2θ = 0.2 rad, radius R = 0.5 m, and thickness h = 0.1 m. The FGM panel is assumed to be with the curved edges fixed and the straight edges simply-supported (FSFS). The panel is under uniaxial compression applied on the curved edges. A 20×20 finite element mesh is used (800 triangular elements). The results for the first three buckling loads parameter ${\bar{N}}_{cr}=N_{cr}R/E_m{h}^2$ are given in Table 5a for the different power law index considered. From this Table 5, a good agreement with the results obtained by Zhao and Liew [22] is observed for the first buckling load, but with discrepancies varying from 4.5% to 9.5 %, for the second and third buckling loads, respectively. This discrepancies can be related to the different models – HSDT and classical finite element used by present model, and FSDT and mesh-free method used by Zhao et al.

Hinged and clamped FGM cylindrical shell panels are now considered to obtain the critical loads. The geometry and mechanical properties are the same of the subsection 5.2. The results for the non-dimensional buckling loads, ${\lambda }_{cr}={N_y}_{cr}R/E_m{H}^2$, obtained using the present model are shown in Table 5b,c, and it is observed that ${\lambda }_{cr}$ decreases when the gradient index p increases. It happens because Ec > EFGM > Em.

5.6 Free Vibrations of Hinged and Clamped FGM Cylindrical Panels

First, for comparison purposes, a clampled FGM cylindrical panel with the following geometry: R/10, a/h=100, E_c=380.0×109 N/m²,${\nu }_c=0.3$, ρ_c=3000 kg/m³, E_m = 70.0×109 N/m², ${\nu }_m=0.3$, ρ_c=2707 kg/m³. The non-dimensional fundamental frequency is given by $\bar{\omega }=\omega a^2\sqrt{{\rho }_m/D_m}$, with $D_m=E_m{H}^3/12(1-{\nu }_m^2)$. The results obtained with the present model are presented and compared in Table 6, with good agreement, with the alternative solutions of Pradyumna and Bandyopadhyay [23] obtained using a higher-order formulation and a C₀ finite element, and the solutions published by Ana et al. [24] obtained using a higher-order shear deformation theory and radial basis functions.

The natural vibrations of hinged and clamped FGM cylindrical shells of section 5.2, with the same geometry and material properties, are now analysed. The first three natural frequencies are obtained for the initial geometry (unloaded panel) considering different gradient index p. Next, the cylindrical panel is loaded by three kinds of loads.

1. Hinged FGM panel loaded:

a) by a uniform pressure on the curved side, from p_y=0 to a final p_y=20.0 ; 40.0 MPa.

b) by a centre point load, from P_c=0 to a final P_c=20.0 ; 40 kN

c) by a uniform pressure on the curved side, from p_y=0 to a final p_y=20.0 ; 40 MPa, combined with a centre point load, from P_c=0 to a final P_c=20.0 ; 40.0 kN.

The geometrically non-linear behaviour is considered, and an incremental-iterative process is used. The results obtained with present model are shown in Table 7a,b. From this Table, it is observed the following:

a) the uniform compressive load increases a little the first two natural frequencies, but decrease a little the third natural frequency.

b) the transverse centre point load decreases more the first three natural frequencies.

c) when both types of loads are applied, the first two frequencies decrease, but less than when only the centre point load is applied.

d) also, it is shown that the frequencies decrease/increase more with the increase of loadings.

2. Clamped FGM panel loaded:

a) by a uniform pressure on the curved side, from p_y=0 to a final p_y=25.0 ; 50.0 MPa

b) by a centre point load, from P_c=0 to a final P_c=20.0; 40.0 kN, or by an external radial pressure, form p_r=0 to a final p_r=2 MPa

c) by a uniform pressure on the curved side, from p_y=0 to a final p_y=25.0 ; 50.0 kPa, combined with a centre point load, from P_c=0 to a final P_c=20.0 ; 40.0 kN, or combined with a external radial pressure, from p_r=0 to a final p_r=1 MPa ; 2 MPa

The corresponding results are shown in Table 8a,b and Table 9a,b. From the Table 8a-b, it is observed the following:

a) the uniform compressive load decreases the first three natural frequencies

b) the transverse centre point load decrease more the first three natural frequencies.

c) when both types of loads are applied, the first three frequencies decrease much more.

d) Also, it is shown that the frequencies decrease more with the increase of loadings.

By comparing the Table 9a,b, it is observed that the external radial pressure decreases more the natural frequencies than the centre point load, as well as their combinations.

Considering the gradient index p=1.0, the variation of the fundamental frequency versus applied type load and boundary conditions, is shown in Figure 10a–c for the case of hinged cylindrical panels, and Figure 11a–c for the case of clamped cylindrical panels. For the case of hinged panel the fundamental frequency is very low when the centre point load is very near the limit load, but increases under axial pressure. For the uniaxial pressure the linear buckling obtained has been 279.6 MPa. However, using the nonlinear analysis the maximum transverse displacement is obtained at load 200 MPa. For the case of the clamped panel, the fundamental frequency first decreases, but after increases in all types of loading.

6. Conclusions

A finite element model for the free vibration analysis of functionally graded material plates/shell panels under mechanical loading is presented. The loading is done by an incremental-iterative process to consider the geometrically nonlinear deformation of the structures.

The finite element model is based on the Reddy’s third-order shear deformation theory applied to a non-conforming triangular flat plate/shell element with 3 nodes, and 8 degrees of freedom per node.

From the present applications it is observed that the FGM plate/shell structures have a smooth variation of stresses through the thickness, and the free vibrations are influenced by the applied loads:

For clamped FGM plates, the uniaxial compressive load leads to a decreasing of the natural frequencies, while the uniform transverse load leads to an increasing of the frequencies. For clamped FGM cylindrical panels, both type of loads lead to a decreasing of the natural frequencies of the free vibration. In contrary, for hinged FGM cylindrical panels, the uniform uniaxial compressive load leads to an increasing of the first two frequencies, while the centre point load decreases the first three frequencies.

From these observations, it can be concluded that loads applied to plate/shell FGM structures can have significant impact on the natural frequencies, depending of the magnitude of the applied loads.

The authors dedicate this paper to the memory of Professor Rolands Rikards (1942- 2022), who sadly passed away on February 2, 2022. He was World Academic and Distinguished Researcher Expert in the fields of Characterization of Material Properties Multidisciplinary Optimization, Solid Mechanics, Finite Element Techniques, Advanced Composite Materials, and Computational Methods. We extend our sorry to his wife and family and Colleagues at Riga Technical University and to his Friends in Latvia.