Investigation of moderately large amplitude effect on nonlinear free vibrations of porous cylinders using multiple scale method

Abstract

This study focuses on the nonlinear free vibrations analysis of porous cylindrical shells for general porosity distribution through thickness. The motion equations which are a system of four coupled nonlinear partial differential equations are extracted in the framework of Mirsky–Herrmann’s shear deformation theory and nonlinear von Karman relations. By employing the multiple-scale method, the system of equations is solved and an analytical solution is proposed for different orders of the equations. By utilizing a parametric study, the effect of different parameters on the nonlinear free vibration behavior of the porous cylinder is discussed. Moreover, the effect of different combinations of boundary conditions and different types of porosity distributions on the linear/nonlinear frequencies is studied. The results are compared with other references in the special cases. The nonlinear frequency behavior of porous cylinders may be hardening or softening and it depends on different conditions, for example, the porosity distribution and boundary conditions. The porosity affects significantly the natural frequency, allowing for substantial weight reduction without compromising frequency in certain cases.

1. Introduction

Porous materials have become a crucial and interesting topic of engineering research due to their unique combination of properties. They offer high strength-to-weight ratios, excellent energy and noise absorption capabilities, and possess a wide range of applications in fields like aerospace, automotive, and biomedical engineering. There are a lot of works in this field for example Wang and Wu (2017) performed a linear free vibration analysis of a functionally graded porous (FGP) cylindrical shell using a sinusoidal shear deformation theory and Rayleigh–Ritz method. Gao et al. (2018) proposed an analytical method for the nonlinear primary resonance analysis of FGP cylindrical shells under uniformly distributed harmonic load. The governing equations were determined using the Donnell shell theory and von Karman relations and by employing the Galerkin method, Duffing-type equations were derived and solved using the multiple-scale method. Li et al. (2019) provided a solution for the linear free vibration analysis of FGP cylindrical shells according to the first-order shear deformation theory (FSDT) and by employing the Fourier series and Rayleigh-Ritz method. Wang et al. (2019) performed the nonlinear vibration analysis of metal-foam cylinders reinforced with graphene platelets. The cylinder was modeled using the improved Donnell nonlinear shell theory and the equations were solved by employing the Galerkin and multiple-scale method for the simply-supported boundary conditions. Daemi and Eipakchi (2019) used the elementary ordinary differential equation theory for the linear free vibrations of FGP cylindrical shells based on the FSDT. Keleshteri and Jelovica (2020) studied the nonlinear free vibrations of FGP cylindrical panels using the FSDT and higher-order shear deformation theory (HSDT) and nonlinear von Karman relations. The nonlinear equations were solved by utilizing the generalized differential quadrature method (DQM). Yang et al. (2021) used the FSDT and von Karman relations to study the nonlinear vibrations of a conical shell truncated by FG graphene platelet-reinforced composite. The natural frequencies were determined using the Galerkin method and the harmonic balance technique. The buckling and nonlinear free vibration analysis of an FGP Euler–Bernoulli micro-beam were presented by Dang and Nguyen (2022) based on the nonlocal strain gradient theory. The solution was proposed using the Galerkin and equivalent linearization methods. Liu et al. (2021) used the improved Donnell nonlinear shell theory and Hamilton’s principle to study the nonlinear vibrations of sandwich FGP cylindrical shells. The nonlinear equations were solved using the Galerkin and multiple-scale methods. Eipakchi and Mahboubi Nasrekani (2021) presented an analytical solution for the nonlinear free vibrations of a sandwich cylindrical shell with a honeycomb core layer based on the FSDT, von Karman relations, and multiple-scale method. Mouthanna et al. (2022); Mouthanna et al. (2023) investigated the nonlinear vibrations of an FGP cylindrical panel with and without eccentrically stiffened using the classical shell theory (CST) and von Karman–Donnell relations. The solutions were determined using the Galerkin method, the airy stress functions, and the fourth-order Runge-Kutta approach. Karimiasl and Alibeigloo (2023) presented an investigation for the nonlinear dynamic behavior of a composite sandwich panel with a piezoelectric layer and a honeycomb core. The nonlinear equations were extracted using the HSDT and von Karman relations and they were solved by employing DQM. Cho (2023) investigated the free vibrations of FGP panels reinforced with graphene platelets numerically using the FSDT and Halpin–Tsai model and by employing the finite element (FE) method. Li et al. (2022) studied the nonlinear free vibrations of FG fiber-reinforced honeycomb sandwich cylindrical shells based on the FSDT, von Karman relations, Halpin–Tsai model, and modified Gibson’s formula. The solution procedure was determined using the Galerkin method, static condensation technique, and multiple-scale method. Ni et al. (2023) employed the Taylor series expansion and DQM to numerically study the nonlinear vibrations of FG graphene nanoplatelet reinforced composite membrane. Pham et al. (2023) proposed a two-node beam element to solve the free vibration problem of an FGP curved nanobeam in a hygro-thermal-magnetic environment based on the Lagrange and Hermite interpolation functions. Shadmani et al. (2023); Shadmani et al. (2024) studied the nonlinear vibrations of temperature-dependent/nondependent truncated conical FG shell based on the FSDT and von Karman relations. The nonlinear equations were solved by employing the Galerkin method and modified Poincare–Lindstedt method. Using the FE method and artificial neural network, a model was proposed by Mahesh (2024) to study the nonlinear free vibrations of FGP plates based on the HSDT and von Karman relations. Ebrahimi et al. (2024) analyzed the nonlinear vibrations of porous truncated conical shells based on the FSDT and von Karman relations. The nonlinear equations were solved using the generalized DQM. The free vibrations of FGP doubly curved panels reinforced with graphene platelets were studied by Zhang et al. (2024) based on the FSDT, Halpin–Tsai model, and the FE method. Amir et al. (2024) studied the free vibrations of FG panels with geometric nonlinearity, porosity, and curvature characteristics by employing the nonlinear FE formulation based on the HSDT and von Karman relations. Wang et al., 2023a analyzed the nonlinear forced vibration characteristics of FG conical shells by employing the FSDT and von Karman relations to extract the governing equations. The equations were solved using the Galerkin and multiple-scale methods. Cong et al. (2022) investigated the nonlinear free vibrations of temperature-dependent FG-CNTRC double-curved shells based on the FSDT and von Karman relations. They used the Airy stress function and Bubnov–Galerkin method to solve the equations. A novel dynamic system for porous cylindrical panels was presented by Li et al. (2023) using the artificial spring technique and based on the FSDT. By employing the DQM, a solution for the free vibrations of the panel with arbitrary boundary conditions was proposed. Wang et al. (2023b) proposed a viscous damping model for the free vibration response analysis based on time-domain damping and eliminating the divergent term for the complex damping model. Kong et al. (2024) presented an innovative negative Poisson’s ratio arc curve resonator structure using a honeycomb core to achieve perfect low-frequency absorption. The structure was modeled using the FE method and validated by experiments. An analysis of the nonlinear dynamic response of a cantilever multilayer microplate with geometric imperfection was conducted by Wu et al. (2024) based on the modified couple stress theory, FSDT, and von Karman relations. The solution was extracted using the Galerkin method. Zhao et al. (2024) employed the nonlinear energy sinks to control the vibration of a single- and double-beam system. Hao et al. (2024) studied the vibrations of rotating porous joined conical–conical shells with elastic supports based on the FSDT and employing the generalized DQM. An investigation of vibrations for variable-walled thickness graphene reinforced porous joined conical–conical panels was conducted by Li et al. (2024) by employing the FSDT and von Karman relations. The governing equations were solved using DQM. Yang et al. (2024) established a new model for FG graphene platelets reinforced composite conical shells under the combined transverse and in-plane forces based on the FSDT and von Karman relations. The results were extracted using the Runge-Kutta method. Jiammeepreecha et al. (2024) presented a comparison between FSDT and HSDT theories for FG spherical and elliptical shells free vibrations under nonlinear temperature distributions by employing the FE method with quadrilateral Lagrangian elements. Xu et al. (2024) used the Gegenbauer–Ritz method to present a solution for the free vibration analysis of rotating FG graphene-reinforced porous cylindrical shells based on the Halpin–Tsai model and the open cell body theory. Aghamaleki et al. (2024) optimized the geometric parameters of FG cylindrical shells resting on an elastic foundation by employing the FSDT and FE method. Hasan and Ali (2024) developed an analytical model for the forced vibrations of composite cylindrical shells with viscous damping and thermal environment effects based on the extended Halpin–Tsai technique, shear deformation theory and von Karman relations. The solution was determined using the Galerkin procedure and the multiple scale method. Wang et al. (2025) utilized the nonlinear energy sink to suppress the vibrations of a variable thickness porous conical shell based on the FSDT and Galerkin method.

This manuscript studies the nonlinear free vibration analysis of FGP cylinders. The governing equations are determined using Mirsky–Herrmann’s shear deformation theory and nonlinear von Karman relations. We employ the multiple-scale analytical method to investigate the effect of material properties and geometric configurations on the nonlinear vibrations of these cylinders. By considering three distinct porosity distributions (symmetric, asymmetric, and uniform) through the thickness, this study offers valuable insights into the influence of porosity variation on the nonlinear free vibrations of FGP cylinders. Although the results are presented for three distinct porosity distributions, the presented solution applies to any type of porosity.

2. Governing equations

Figure 1 shows an FGP cylinder with length L, thickness h, and its coordinate system. The origin of the cylindrical coordinate system has been located on the mid-plane and R _m is the radius of the mid-plane. The location of each point of the FGP cylinder is defined using three parameters r, x, and θ, where r = R _m +z and z are measured from the mid-plane.

Figure 1.

FGP cylinder and coordinate system.

By assuming axisymmetric condition and according to Mirsky–Herrmann’s shear deformation theory, the displacement field has been considered as follows (Mirsky and Herrmann, 2021):

$\begin{array}{l}{U}_x\left(z,x,t\right)=u\left(x,t\right)+z.\phi \left(x,t\right)\\ U_z\left(z,x,t\right)=w\left(x,t\right)+z.\psi \left(x,t\right);U_{\theta }\left(z,x,t\right)=0\end{array}$

(1)

where U _x , U _z , and U _θ are displacements in x, z, and θ directions, respectively, and u, φ, w, and ψ are unknown functions. The strain displacement relations for moderately large deformations are considered according to the von Karman relations as follows (Sadd, 2009):

$\begin{array}{l}{\epsilon }_x=\frac{\partial U_x}{\partial x}+\frac{1}{2}{\left(\frac{\partial U_z}{\partial x}\right)}^2=\frac{\partial u}{\partial x}+z\frac{\partial \phi }{\partial x}+\frac{z^2}{2}{\left(\frac{\partial \psi }{\partial x}\right)}^2+\frac{\partial w}{\partial x}\frac{\partial \psi }{\partial x}+\frac{1}{2}{\left(\frac{\partial w}{\partial x}\right)}^2\\ {\epsilon }_{\theta }=\frac{U_z}{r}=\frac{w+z\psi }{R_m+z};\hspace{0.5em}{\epsilon }_z=\frac{\partial U_z}{\partial z}+\frac{1}{2}{\left(\frac{\partial U_z}{\partial z}\right)}^2=\psi ·\left(1+\frac{\psi }{2}\right);\hspace{0.5em}{\gamma }_{x\theta }={\gamma }_{z\theta }=0\\ {\gamma }_{xz}=\frac{\partial U_x}{\partial z}+\frac{\partial U_z}{\partial x}+\frac{\partial U_z}{\partial x}·\frac{\partial U_z}{\partial z}=\phi +z\frac{\partial \psi }{\partial x}.\left(1+\psi \right)+\frac{\partial w}{\partial x}.\left(1+\psi \right)\end{array}$

(2)

The stress-strain relations for an FGP cylinder are defined based on the Biot theory as follows (Biot, 1962):

$\begin{array}{l}{\sigma }_x=\left(K-2G/3\right)\hspace{0.17em}\left({\epsilon }_x+{\epsilon }_{\theta }+{\epsilon }_z\right)+2G\hspace{0.17em}{\epsilon }_x;\hspace{1em}{\sigma }_{\theta }=\left(K-2G/3\right)\hspace{0.17em}\left({\epsilon }_x+{\epsilon }_{\theta }+{\epsilon }_z\right)+2G\hspace{0.17em}{\epsilon }_{\theta }\\ {\sigma }_z=\left(K-2G/3\right)\hspace{0.17em}\left({\epsilon }_x+{\epsilon }_{\theta }+{\epsilon }_z\right)+2G\hspace{0.17em}{\epsilon }_z;\hspace{1em}{\tau }_{xz}=G\hspace{0.17em}{\gamma }_{xz};\hspace{1em}{\tau }_{x\theta }={\tau }_{z\theta }=0\end{array}$

(3)

where K and G are bulk and shear modulus, respectively, and they are functions of z and considered as K = K₀ λ(z) and G = G₀ λ(z) in which λ(z) is a dimensionless function of z and depends on the porosity distributions (symmetric, asymmetric, and uniform). K₀ and G₀ are the maximum/solid bulk and shear modulus without porosity. Poisson’s ratio is considered constant and to be unaffected by the porosity. It should be mentioned that stress-strain relations are derived from the assumption that the fluid in the porous material is air (with low pressure) and the cavity pressure is zero in this study. The strain energy and its variation are determined as follows:

$U=\underset{V}{\iiint }{U}^{*}dV,\hspace{1em}dV=R_m\left(1+\frac{z}{R_m}\right)\hspace{0.5em}d\theta \hspace{0.5em}dz\hspace{0.5em}dx,\hspace{1em}0\le \theta \le 2\pi ,-\frac{h}{2}\le z\le \frac{h}{2},\hspace{1em}0\le x\le L$

(4a)

$U^{*}=\frac{1}{2}\left({\sigma }_x\hspace{0.17em}{\epsilon }_x+{\sigma }_{\theta }\hspace{0.17em}{\epsilon }_{\theta }+{\sigma }_z\hspace{0.17em}{\epsilon }_z+{\tau }_{xz}\hspace{0.17em}{\gamma }_{xz}\right)$

(4b)

$\frac{\delta U}{2\pi }=R_m\hspace{0.5em}\underset{0}{\overset{L}{{\displaystyle \int }}}\underset{-h/2}{\overset{h/2}{{\displaystyle \int }}}\left({\sigma }_x\hspace{0.5em}\delta {\epsilon }_x+{\sigma }_{\theta }\hspace{0.5em}\delta {\epsilon }_{\theta }+{\sigma }_z\hspace{0.5em}\delta {\epsilon }_z+{\tau }_{xz}\hspace{0.5em}\delta {\gamma }_{xz}\right)\left(1+\frac{z}{R_m}\right)\hspace{0.5em}dz\hspace{0.5em}dx$

(4c)

where h and L are the thickness and length of the cylinder. Stress resultants are defined as follows:

$\begin{array}{l}{\left\{\begin{array}{lll}{N}_x& M_x& Q_x\end{array}\right\}}^T={\int }_{-\frac{h}{2}}^{\frac{h}{2}}{\sigma }_x{\left\{\begin{array}{lll}1& z& z^2\end{array}\right\}}^T\left(1+\frac{z}{R_m}\right)\hspace{0.17em}dz;\hspace{0.5em}{N}_z={\int }_{-h/2}^{h/2}{\sigma }_z\left(1+\frac{z}{R_m}\right)\hspace{0.17em}dz\\ {\left\{\begin{array}{ll}{N}_{\theta }& M_{\theta }\end{array}\right\}}^T={\int }_{-h/2}^{h/2}{\sigma }_{\theta }{\left\{\begin{array}{ll}1& z\end{array}\right\}}^{T}dz;\hspace{1em}{\left\{\begin{array}{ll}{N}_{xz}& M_{xz}\end{array}\right\}}^T={\kappa }_s{\int }_{-h/2}^{h/2}{\tau }_{xz}\left\{\begin{array}{ll}1& z\end{array}\right\}\left(1+\frac{z}{R_m}\right)\hspace{0.17em}dz\end{array}$

(5)

where κ _s is the shear correction factor which is considered 5/6 in this study (Vlachoutsis, 1992). By substituting equations (1)–(4) into (5), the stress resultants can be extracted as functions of the displacement field. For example, N _x is as follows:

$\begin{array}{l}{N}_x={\int }_{-\frac{h}{2}}^{\frac{h}{2}}{\sigma }_x\left(1+\frac{z}{R_m}\right)\hspace{0.17em}dz=a_1{J}_{01}\frac{\partial u}{\partial x}+a_1{J}_{11}\frac{\partial \phi }{\partial x}+\frac{1}{2}{a}_1{J}_{01}{\left(\frac{\partial w}{\partial x}\right)}^2+\frac{b_2}{R_m}{J}_{00}w+a_1{J}_{11}\frac{\partial w}{\partial x}\frac{\partial \psi }{\partial x}\\ +\frac{1}{2}{a}_1{J}_{21}{\left(\frac{\partial \psi }{\partial x}\right)}^2+\frac{1}{2}{b}_2{J}_{01}{\psi }^2+b_2\left(J_{01}+\frac{J_{10}}{R_m}\right)\psi \end{array}$

(6a)

$a_1=\frac{4}{3}{G}_0+K_0;\hspace{1em}{b}_1=\frac{2}{3}{G}_0+K_0;\hspace{1em}{b}_2=-\frac{2}{3}{G}_0+K_0;\hspace{1em}{J}_{i,j}={\int }_{-\frac{h}{2}}^{\frac{h}{2}}{z}^i\lambda \left(z\right).{\left(1+\frac{z}{R_m}\right)}^{j}dz$

(6b)

The kinetic energy for the FGP cylinder is defined as follows:

$T=\underset{V}{\iiint }{T}^{*}dV,dV=R_m\hspace{0.17em}\left(1+\frac{z}{R_m}\right)\hspace{0.5em}d\theta \hspace{0.5em}dz\hspace{0.5em}dx,0\le \theta \le 2\pi ,-\frac{h}{2}\le z\le \frac{h}{2},\hspace{1em}0\le x\le L$

(7a)

$T^{*}=\frac{\rho \left(z\right)}{2}\hspace{0.5em}\left({\left(\frac{\partial U_x}{\partial t}\right)}^2+{\left(\frac{\partial U_z}{\partial t}\right)}^2\right)$

(7b)

$\begin{array}{l}\frac{\delta T}{2\pi }=-{\rho }_s{\int }_0^L[\left({\mathrm{{\rm I}}}_0\frac{{\partial }^{2}u}{\partial t^2}+{\mathrm{{\rm I}}}_1\frac{{\partial }^2\phi }{\partial t^2}\right)\delta \hspace{0.17em}u+\left({\mathrm{{\rm I}}}_1\frac{{\partial }^{2}u}{\partial t^2}+{\mathrm{{\rm I}}}_2\frac{{\partial }^2\phi }{\partial t^2}\right)\delta \hspace{0.17em}\phi +\left({\mathrm{{\rm I}}}_0\frac{{\partial }^{2}w}{\partial t^2}+{\mathrm{{\rm I}}}_1\frac{{\partial }^2\psi }{\partial t^2}\right)\delta \hspace{0.17em}w\\ +\left({\mathrm{{\rm I}}}_1\frac{{\partial }^{2}w}{\partial t^2}+{\mathrm{{\rm I}}}_2\frac{{\partial }^2\psi }{\partial t^2}\right)\delta \hspace{0.17em}\psi ]dx\\ {\mathrm{{\rm I}}}_n={\int }_{-\frac{h}{2}}^{\frac{h}{2}}{z}^n\hspace{0.5em}\mathrm{\Gamma }\left(z\right)\left(1+\frac{z}{R_m}\right)dz,\hspace{1em}\rho \left(z\right)={\rho }_s\mathrm{\Gamma }\left(z\right)\end{array}$

(7c)

where ρ(z) is the density of the FGP cylinder and Г(z) is the dimensionless volume function and shows the density variation through the thickness. The dimensionless volume function Г(z) is defined according to the porosity distributions (symmetric, asymmetric, and uniform) which is discussed later. ρ _s is the density of the solid cell in the absence of porosity. Employing the Hamilton principle which states $\delta {\int }_{t_1}^{t_2}\left(T-U\right)dt=0$ the governing equations and boundary conditions are derived as follows (Rao, 2007):

$E{q}_1:{N_x}^{\prime }-{\rho }_s\hspace{0.17em}\left({\mathrm{{\rm I}}}_0\frac{{\partial }^{2}u}{\partial t^2}+{\mathrm{{\rm I}}}_1\frac{{\partial }^2\phi }{\partial t^2}\right)=0;$

(8a)

$E{q}_2:{M_x}^{\prime }-N_{xz}-{\rho }_s\hspace{0.17em}\left({\mathrm{{\rm I}}}_2\frac{{\partial }^2\phi }{\partial t^2}+{\mathrm{{\rm I}}}_1\frac{{\partial }^{2}u}{\partial t^2}\right)=0$

(8b)

$E{q}_3:{\left(N_x\hspace{0.5em}{w}^{\prime }+M_x\hspace{0.5em}{\psi }^{\prime }+N_{xz}\left(\psi +1\right)\right)}^{\prime }-\frac{N_{\theta }}{R}-{\rho }_s\left({\mathrm{{\rm I}}}_0\frac{{\partial }^{2}w}{\partial t^2}+{\mathrm{{\rm I}}}_1\frac{{\partial }^2\psi }{\partial t^2}\right)=0$

(8c)

$E{q}_4:{\left(M_x{w}^{\prime }+Q_x{\psi }^{\prime }+M_{xz}\left(\psi +1\right)\right)}^{\prime }-{\psi }^{\prime }{M}_{xz}-N_z\left(\psi +1\right)-w^{\prime }{N}_{xz}-\frac{M_{\theta }}{R}-{\rho }_s\left({\mathrm{{\rm I}}}_1\frac{{\partial }^{2}w}{\partial t^2}+{\mathrm{{\rm I}}}_2\frac{{\partial }^2\psi }{\partial t^2}\right)=0$

(8d)

$\begin{array}{l}B.C:{\left(N_x\delta \hspace{0.17em}u+M_x\delta \hspace{0.17em}\phi +H_1\delta \hspace{0.17em}w+H_2\delta \hspace{0.17em}\psi \right)|}_0^L=0;\hspace{1em}{\mathrm{A}}^{\prime }=\frac{\partial \mathrm{A}}{\partial x}\\ H_1=N_x\hspace{0.5em}{w}^{\prime }+M_x\hspace{0.5em}{\psi }^{\prime }+N_{xz}\left(\psi +1\right);\hspace{1em}{H}_2=M_x\hspace{0.5em}{w}^{\prime }+Q_x\hspace{0.5em}{\psi }^{\prime }+M_{xz}\left(\psi +1\right)\end{array}$

(8e)

In this study, the boundary conditions for special cases are considered as follows:

$\text{Clamped}\hspace{0.17em}\text{edge}:u=0,w=0,\phi =0,\psi =0$

(9a)

$\text{Simply}\hspace{0.17em}\text{supported}\hspace{0.17em}\text{edge}\hspace{0.17em}\left(\text{case}\hspace{0.17em}1\right):\hspace{1em}u=0,w=0,M_x=0,H_2=0$

(9b)

$\text{Simply}\hspace{0.17em}\text{supported}\hspace{0.17em}\text{edge}\hspace{0.17em}\left(\text{case}\hspace{0.17em}2\right):\hspace{1em}{N}_x=0,w=0,M_x=0,H_2=0$

(9c)

$\text{Free}\hspace{0.17em}\text{edge}:\hspace{1em}{N}_x=0,M_x=0,H_1=0,H_2=0$

(9d)

where H₁ and H₂ are defined in equation (8d). For the free vibration analysis and simply supported boundary conditions, in this study, case 2 is considered. By substituting equations (1)–(3), and (5) into equation (7), the equations of motion and boundary conditions are determined as functions of displacement field parameters as follows:

$\begin{array}{l}E{q}_1:{\mathrm{\Pi }}_1\left(u,\phi ,w,\psi ,x,t,\partial /\partial x,{\partial }^2/\partial x^2\right)-{\rho }_s\hspace{0.17em}\left({\mathrm{{\rm I}}}_0\frac{{\partial }^{2}u}{\partial t^2}+{\mathrm{{\rm I}}}_1\frac{{\partial }^2\phi }{\partial t^2}\right)=0\\ {\mathrm{\Pi }}_1=a_1\left(J_{01}\frac{{\partial }^{2}u}{\partial x^2}+J_{11}\frac{{\partial }^2\phi }{\partial x^2}\right)+a_1\left(J_{01}\frac{\partial w}{\partial x}+J_{11}\frac{\partial \psi }{\partial x}\right)\frac{{\partial }^{2}w}{\partial x^2}+\frac{b_2{J}_{00}}{R}\frac{\partial w}{\partial x}+b_2\left(J_{01}\left(\psi +1\right)+\frac{J_{10}}{R}\right)\frac{\partial \psi }{\partial x}\\ +a_1\left(J_{11}\frac{\partial w}{\partial x}+J_{21}\frac{\partial \psi }{\partial x}\right)\frac{{\partial }^2\psi }{\partial x^2}\end{array}$

(10a)

$\begin{array}{l}E{q}_2:{\mathrm{\Pi }}_2\left(u,\phi ,w,\psi ,x,t,\partial /\partial x,{\partial }^2/\partial x^2\right)-{\rho }_s\hspace{0.17em}\left({\mathrm{{\rm I}}}_2\frac{{\partial }^2\phi }{\partial t^2}+{\mathrm{{\rm I}}}_1\frac{{\partial }^{2}u}{\partial t^2}\right)=0\\ {\mathrm{\Pi }}_2=a_1\hspace{0.5em}\left(\hspace{0.5em}{J}_{11}\hspace{0.5em}\frac{{\partial }^{2}u}{\partial x^2}+J_{21}\hspace{0.5em}\frac{{\partial }^2\phi }{\partial x^2}\right)+a_1\left(J_{11}\hspace{0.5em}\frac{\partial w}{\partial x}+J_{21}\hspace{0.5em}\frac{\partial \psi }{\partial x}\right)\frac{{\partial }^{2}w}{\partial x^2}+a_1\hspace{0.5em}\left(J_{21}\frac{\partial w}{\partial x}+J_{31}\frac{\partial \psi }{\partial x}\right)\hspace{0.5em}\frac{{\partial }^2\psi }{\partial x^2}\\ +\left(-{\kappa }_s{J}_{01}{G}_0\left(\psi +1\right)+\frac{b_2{J}_{10}}{R}\right)\frac{\partial w}{\partial x}+\left(J_{11}\left(-{\kappa }_s{G}_0+b_2\right)\left(\hspace{0.17em}\psi +1\right)+\frac{b_2{J}_{20}}{R}\right)\frac{\partial \psi }{\partial x}-{\kappa }_s{J}_{01}{G}_0\phi \end{array}$

(10b)

$E{q}_3:{\mathrm{\Pi }}_3\left(u,\phi ,w,\psi ,x,t,\partial /\partial x,{\partial }^2/\partial x^2\right)-{\rho }_s\left({\mathrm{{\rm I}}}_0\frac{{\partial }^{2}w}{\partial t^2}+{\mathrm{{\rm I}}}_1\frac{{\partial }^2\psi }{\partial t^2}\right)=0$

(10c)

$E{q}_4:{\mathrm{\Pi }}_4\left(u,\phi ,w,\psi ,x,t,\partial /\partial x,{\partial }^2/\partial x^2\right)-{\rho }_s\left({\mathrm{{\rm I}}}_1\frac{{\partial }^{2}w}{\partial t^2}+{\mathrm{{\rm I}}}_2\frac{{\partial }^2\psi }{\partial t^2}\right)=0$

(10d)

where Π₃ and Π₄ are differential operators similar to II_{1 and} Π₂ which were not reported here due to their long forms.

3. Analytical solution

The equations of motion for an FGP cylinder (equation (10)) are a system of four coupled nonlinear partial differential equations. In this study, the multiple-scale method has been used to solve these equations. For this aim, the following dimensionless parameters are used to determine the dimensionless form of equations:

$x^{*}=\frac{x}{L},t^{*}=\frac{t}{t_0},R^{*}=\frac{R}{R_0},h^{*}=\frac{h}{h_0},u^{*}=\frac{u\hspace{0.17em}\left(x,t\right)}{h_0},w^{*}=\frac{w\hspace{0.17em}\left(x,t\right)}{h_0}$

(11a)

where h₀, R₀, and t₀ are thickness, radius, and time characteristics, respectively, and in this study h₀ = h, R₀ = R _m . By substituting equation (11a) into equation (10), the following dimensionless quantities appeared too:

${\mathrm{{\rm I}}}_n^{*}=\frac{{\mathrm{{\rm I}}}_n}{{h_0}^{n+1}},J_{n,m}^{*}=\frac{J_{n,m}}{{h_0}^{n+1}},\beta =\frac{G_0}{K_0},\epsilon =\frac{h_0}{L},e=\frac{{\rho }_s{h_0}^2}{K_0{t_0}^2},z_2=\frac{R_0}{h_0}$

(11b)

We defined t₀ = L/c, where $c=\sqrt{K_0/{\rho }_s}$ . ε is the small quantity which is considered as the perturbation parameter. By increasing the length or decreasing the thickness, this parameter decreases. It should be mentioned that φ and ψ are dimensionless parameters. To solve the equations a new variable X=(x*−1)/ε is defined and applied to the equations. To solve the equations and investigate the mutual effect of amplitude and frequency, the multiple-scale method is a suitable approach (Nayfeh, 1993). To apply this method, new time scales including fast and slow time scales T₀ = t*, T₁ = εt*, T₂ = ε²t*, … are defined. By employing the new time scales and according to the chain rule, the time derivatives in the equations are substituted by the following derivative (Nayfeh, 1993):

$\frac{{\partial }^2}{\partial t^2}=D_0^2+2\epsilon \hspace{0.17em}{D}_0{D}_1+{\epsilon }^2\left(2D_0{D}_2+D_1^2\right)+2{\epsilon }^3{D}_1{D}_2+...;\hspace{1em}{D}_n=\frac{\partial }{\partial T_n},n=0,1,2,...$

(12)

By applying the new dimensionless parameters, variables, and new time scales, the dimensionless form of equations is determined. For example, the first equation is as follows:

$\begin{array}{l}E{q}_1:{\mathrm{\Pi }}_1^{*}\left(u^{*},\phi ,w^{*},\psi ,X,T_0,T_1,T_2,\epsilon \right)=0;\hspace{1em}{A}_{11}=I_0^{*}u+I_1^{*}\phi \\ {\mathrm{\Pi }}_1^{*}=\left(3-2\beta \right)\left(z_2{R}^{*}{J}_{01}^{*}\psi \frac{\partial \psi }{\partial X}+J_{00}^{*}\frac{\partial w}{\partial X}\right)+\left(z_2{R}^{*}\left(3+4\beta \right)\left(J_{01}^{*}\frac{{\partial }^{2}w}{\partial X^2}+J_{11}^{*}\frac{{\partial }^2\psi }{\partial X^2}\right)\right)\frac{\partial w}{\partial X}+\left(3-2\beta \right)(J_{10}^{*}\\ +z_2{R}^{*}{J}_{01}^{*}))\frac{\partial \psi }{\partial X}+(z_2{R}^{*}\left(3+4\beta \right)\left(J_{11}^{*}\frac{{\partial }^{2}w}{\partial X^2}+J_{21}^{*}\frac{{\partial }^2\psi }{\partial X^2}\right)+z_2{R}^{*}\left(3+4\beta \right)\left(J_{01}^{*}\frac{{\partial }^{2}u}{\partial X^2}+J_{11}^{*}\frac{{\partial }^2\phi }{\partial X^2}\right)\\ -3z_2{R}^{*}e\left(\frac{{\partial }^2{A}_{11}}{\partial {T_0}^2}+2\epsilon \frac{{\partial }^2{A}_{11}}{\partial T_0\partial T_1}+{\epsilon }^2\frac{{\partial }^2{A}_{11}}{\partial {T_1}^2}+2{\epsilon }^2\frac{{\partial }^2{A}_{11}}{\partial T_0\partial T_2}\right)\end{array}$

(13)

The response is considered as the following expansion:

$\begin{array}{l}{u}^{*}\hspace{0.17em}\left(X,T_0,T_1,T_2\right)=\epsilon .{\displaystyle \sum }_{i=0}^2{\epsilon }^i{u}_i\hspace{0.17em}\left(X,T_0,T_1,T_2\right);\hspace{0.5em}{\phi }^{*}\hspace{0.17em}\left(X,T_0,T_1,T_2\right)=\epsilon .{\displaystyle \sum }_{i=0}^2{\epsilon }^i{\phi }_i\hspace{0.17em}\left(X,T_0,T_1,T_2\right)\\ w^{*}\hspace{0.17em}\left(X,T_0,T_1,T_2\right)=\epsilon .{\displaystyle \sum }_{i=0}^2{\epsilon }^i{w}_i\hspace{0.17em}\left(X,T_0,T_1,T_2\right);\hspace{0.5em}{\psi }^{*}\hspace{0.17em}\left(X,T_0,T_1,T_2\right)=\epsilon .{\displaystyle \sum }_{i=0}^2{\epsilon }^i{\psi }_i\hspace{0.17em}\left(X,T_0,T_1,T_2\right)\end{array}$

(14)

By substituting equation (14) into dimensionless equations and separating terms with the same order of ε, equations of order zero, and one, etc. are determined. The equations of order-zero are as follows:

$\begin{array}{l}{\left[B_2\right]}_{4\times 4}\frac{{\partial }^2}{\partial X^2}\left\{y_0\right\}+{\left[B_1\right]}_{4\times 4}\frac{\partial }{\partial X}\left\{y_0\right\}+{\left[B_0\right]}_{4\times 4}\left\{y_0\right\}+{\left[B_3\right]}_{4\times 4}\frac{{\partial }^2}{\partial {T_0}^2}\left\{y_0\right\}={\left\{0\right\}}_{4\ast 1}\\ \left\{y_0\right\}={\left\{\begin{array}{llll}{u}_0& {\phi }_0& w_0& {\psi }_0\end{array}\right\}}^T\end{array}$

(15)

where matrices of coefficients have been reported in the Appendix. The solution of equation (15) which is a system of four coupled linear partial differential equations with constant coefficients is considered as follows:

$\left\{y_0\right\}={\left\{V\left(X,T_1,T_2\right)\right\}}_{4\times 1}\hspace{0.5em}{e}^{\hspace{0.17em}i\omega T_0}$

(16)

where ω and {V(X,T₁,T₂)} are the linear dimensionless natural frequency and corresponding mode shapes, respectively. If one multiplies ω by t ₀ (i.e., ω.t ₀ ) the frequency in terms of rad/s is determined. By substituting equation (16) into equation (15), a system of four coupled partial differential equations is extracted as follows:

$\left[B_2\right]\frac{{\partial }^2}{\partial X^2}\left\{V\right\}+\left[B_1\right]\frac{\partial }{\partial X}\left\{V\right\}+\left[B_0\right]\left\{V\right\}-i.{\omega }^2\left[B_3\right]\left\{V\right\}={\left\{0\right\}}_{4\ast 1};\hspace{1em}\left\{V\right\}={\left\{V\left(X,T_1,T_2\right)\right\}}_{4\times 1}$

(17)

The solution of equation (17) is assumed as follows:

$\left\{V\right\}={\left\{A\left(T_1,T_2\right)\right\}}_{4\times 1}\hspace{0.5em}{e}^{\alpha X}$

(18)

where α and {A(T₁,T₂)} are eigenvalues and corresponding eigenvectors, respectively. By substituting equation (18) into equation (17), a system of four coupled algebraic equations is determined as the following:

$\left[B_T\right]\left\{A\right\}={\left\{0\right\}}_{4\times 1};\hspace{1em}\left\{A\right\}={\left\{A\left(T_1,T_2\right)\right\}}_{4\times 1};\hspace{1em}{\left[B_T\right]}_{4\times 4}=\left[B_2\right]{\alpha }^2+\left[B_1\right]\alpha +\left[B_0\right]-i.{\omega }^2\left[B_3\right]$

(19)

For a nontrivial solution, the determinant of the coefficient matrix must be zero (i.e., det[B _T ] = 0), which results in an algebraic equation from order-eight that is called the dispersion equation, and by solving it, eight eigenvalues α _i , i = 1..8 and eight corresponding eigenvectors {A} _i , i = 1..8 are determined. The general solution is considered as follows:

$\left\{V\right\}={\displaystyle \sum }_{j=1}^8{C}_j{\left\{A\right\}}_j\hspace{0.5em}{e}^{{\alpha }_{j}X};\hspace{1em}{C}_j=C_j\left(T_1,T_2\right)$

(20)

where C _j are constant coefficients which are determined by the boundary conditions. By applying the boundary conditions at two edges, a system of eight coupled algebraic equations is obtained as [M]_8*8{C}_8*1 = {0}_8*1 where {C}_8*1 = {C₁ C₂ C₃ C₄ C₅ C₆ C₇ C₈} ^T and by solving this eigenvalue problem, it is possible to determine all constants as a function of one of the constants (i.e., C _j = c _j *C₁, j = 2,…,8 and c _j are known values). To find the nontrivial solution, det[M] has been equated to zero which results in a complicated algebraic equation with respect to ω. We used the bisection method to solve this equation. Finally, the general form of the solution of order-zero is determined as follows:

$\left\{y_0\right\}=C_1\left(T_1,T_2\right){\displaystyle \sum }_{j=1}^8{\left\{A^{\prime }\right\}}_j\hspace{0.5em}{e}^{{\alpha }_{j}X}{e}^{\hspace{0.17em}i\omega T_0}+C.C.$

(21)

where C.C. stands for the complex conjugate of terms and {A´} _j are the normalized eigenvectors that contain the known values c _j . Now, the solution of order one will be explained. Equations of order one are as follows:

$\begin{array}{l}{\left[B_2\right]}_{4\times 4}\frac{{\partial }^2}{\partial X^2}\left\{y_1\right\}+{\left[B_1\right]}_{4\times 4}\frac{\partial }{\partial X}\left\{y_1\right\}+{\left[B_0\right]}_{4\times 4}\left\{y_1\right\}+{\left[B_3\right]}_{4\times 4}\frac{{\partial }^2}{\partial {T_0}^2}\left\{y_1\right\}=-\left\{F_1\right\}\\ \left\{y_1\right\}={\left\{\begin{array}{llll}{u}_1& {\phi }_1& w_1& {\psi }_1\end{array}\right\}}^T;\hspace{1em}\left\{F_1\right\}=\left\{\begin{array}{llll}{f}_{11}& f_{12}& f_{13}& f_{14}\end{array}\right\}\end{array}$

(22)

The homogeneous part of equations order-one (equation (22)) is similar to equations order zero (equation (15)) and the nonhomogeneous terms are functions of the solution of order-zero (i.e., u₀, φ₀, w₀, and ψ₀) which have been reported in the Appendix. Equation (22) is a system of four coupled linear nonhomogeneous partial differential equations with constant coefficients. In the multiple-scale method, for the higher order equations, only the particular solution is considered and the total solution of the equations is the summation of the solution of order-zero, and the particular solutions of higher orders (i.e., orders-one and two) (Nayfeh, 1993). By substituting the solution of order-zero (equation (21)) into the equations of order-one (equation (22)), the nonhomogeneous terms are determined as follows:

$\left\{F_1\right\}={\left\{F_{11}\right\}}_{4\times 1}\hspace{0.17em}\exp \left(i\omega T_0\right)+{\left\{F_{12}\right\}}_{4\times 1}\hspace{0.17em}\exp \left(2i\omega T_0\right)+C.C$

(23)

where coefficient vectors {F₁₁} and {F₁₂} are known functions of X, T₁, and T₂. As it is seen, the nonhomogeneous terms have two parts: secular terms including exp(iωT₀) and exp(−iωT₀) which should be vanished, and nonsecular terms including exp(2iωT₀) and exp(−2iωT₀) which its corresponding solution should be determined. The solution for the nonsecular terms is considered as follows:

$\left\{y_{p1}\right\}={\left\{A_1\right\}}_{4\times 1}\hspace{0.17em}\exp \left(2i\omega T_0\right)\to \left[B_2\right]\frac{{\partial }^2}{\partial X^2}\left\{A_1\right\}+\left[B_1\right]\frac{\partial }{\partial X}\left\{A_1\right\}+\left(\left[B_0\right]-4{\omega }^2\left[B_3\right]\right)\left\{A_1\right\}=\left\{F_{12}\right\}$

(24a)

where {A₁} is an unknown vector which is a function of X, T₁, and T₂. The general form of {A₁} depends on the form of nonhomogeneous terms and vector {F₁₂}. In this work, the nonhomogeneous terms (i.e., {F₁₂}) include exponential (e.g., exp (β _i X) where β _i are eigenvalues or a combination of the eigenvalues of the homogenous equations) and nonexponential components. Consequently, the general form of the nonhomogeneous solution is considered as {A₁} = {K₀}+{K _i }exp(β _i X), and by substituting this solution into the nonhomogeneous equations, the unknown coefficient vectors {K₀} and {K _i } are calculated. For the secular terms, it is required to apply the solvability conditions (Nayfeh, 1993). For this purpose, the solution for these terms is assumed as the following:

$\left\{y_{p2}\right\}=\left\{A_2\right\}\exp \left(i\omega T_0\right)\to EQ=\left[B_2\right]\frac{{\partial }^2\left\{A_2\right\}}{\partial X^2}+\left[B_1\right]\frac{\partial \left\{A_2\right\}}{\partial X}+\left(\left[B_0\right]-{\omega }^2\left[B_3\right]\right)\left\{A_2\right\}-\left\{F_{12}\right\}=\left\{0\right\}$

(24b)

Equation (24b) includes four coupled equations. The solvability condition is considered as equation (25):

$\underset{0}{\overset{L_0}{{\displaystyle \int }}}\left(EQ\left[1\right].{\varphi }_1+EQ\left[2\right].{\varphi }_2+EQ\left[3\right].{\varphi }_3+EQ\left[4\right].{\varphi }_4\right)dX=0;\hspace{0.5em}{L}_0=L/h_0$

(25)

where ϕ₁, ϕ₂, ϕ₃, and ϕ₄ are adjoint functions and they are determined from the homogeneous form of equations and homogeneous boundary conditions. The procedure of determination of the adjoint functions was explained in detail by Eipakchi and Mahboubi Nasrekani (Eipakchi and Mahboubi Nasrekani, 2021). By employing equation (25), the following equation is extracted:

$\frac{\partial C_1\left(T_1,T_2\right)}{\partial T_1}=0\to C_1\left(T_1,T_2\right)=C_1\left(T_2\right)$

(26)

where C₁ is the coefficient of the solution for order-zero in equation (21). Finally, the solution of order-one which includes only the particular solution due to the nonsecular terms is determined in the following general form:

$\left\{y_1\right\}=\left\{A_0\right\}\exp \left(2i\omega T_0\right)+\left\{\overline{A_0}\right\}\exp \left(-2i\omega T_0\right)$

(27)

The bar sign shows the complex conjugate of the quantity. In the final step, it is required to solve the equations of order two. These equations are as follows:

$\begin{array}{l}{\left[B_2\right]}_{4\times 4}\frac{{\partial }^2}{\partial X^2}\left\{y_2\right\}+{\left[B_1\right]}_{4\times 4}\frac{\partial }{\partial X}\left\{y_2\right\}+{\left[B_0\right]}_{4\times 4}\left\{y_2\right\}+{\left[B_3\right]}_{4\times 4}\frac{{\partial }^2}{\partial {T_0}^2}\left\{y_2\right\}=-\left\{F_2\right\}\\ \left\{y_2\right\}={\left\{\begin{array}{llll}{u}_2& {\phi }_2& w_2& {\psi }_2\end{array}\right\}}^T;\hspace{1em}\left\{F_2\right\}=\left\{\begin{array}{llll}{f}_{21}& f_{22}& f_{23}& f_{24}\end{array}\right\}\end{array}$

(28)

f₂₁, f₂₂, f₂₃, and f₂₄ are known functions of solutions of the previous orders (i.e., u₀, φ₀, w₀, ψ₀, u₁, φ₁, w₁, and ψ₁) and due to their long-form, they are not reported here. By substituting the solution of the previous orders (i.e., equations (21) and 27) into equation (28), the nonhomogeneous part for order two is determined as equation (29):

$\left\{F_2\right\}=\left\{F_{21}\right\}\exp \left(i\omega T_0\right)+\left\{F_{22}\right\}\exp \left(3i\omega T_0\right)+C.C$

(29)

where {F₂₁} and {F₂₂} are known functions of X, T₁, and T₂. Similar to the equations of order one, the nonhomogeneous terms of order two include nonsecular and secular terms. For the secular terms, by applying the solvability conditions, the following equation is derived:

$\frac{d}{d{T}_2}{C}_1\left(T_2\right)+\left({\kappa }_R-{\kappa }_{I}i\right)\hspace{0.5em}\overline{C_1}\left(T_2\right)·C_1{\left(T_2\right)}^2=0$

(30)

where κ _R and κ _I are known real coefficients. To solve equation (30), the solution is assumed as equation (31):

$C_1\left(T_2\right)=\frac{1}{2}a\left(T_2\right).e^{i\beta \left(T_2\right)}\to \hspace{1em}\overline{C_1}\left(T_2\right)=\frac{1}{2}a\left(T_2\right).e^{-i\beta \left(T_2\right)}$

(31)

where a(T₂) and β(T₂) are real functions and they stand for the amplitude and phase, respectively. By substituting equation (31) into, equation (30), and considering the real and imaginary parts, the following equations are obtained:

$\frac{d}{d{T}_2}a\left(T_2\right)=-\frac{1}{4}a{\left(T_2\right)}^3{\kappa }_R;\hspace{1em}\frac{d}{d{T}_2}\beta \left(T_2\right)=\frac{1}{4}a{\left(T_2\right)}^2{\kappa }_I$

(32)

By solving equation (32), and applying the initial conditions the amplitude and phase are obtained as follows:

$I.C:a\left(0\right)=a_0,\beta \left(0\right)={\beta }_0\hspace{0.17em}\to a\left(T_2\right)=\pm \frac{a_0}{\sqrt{1+\frac{{\kappa }_R}{2}{T}_2{a_0}^2}},\beta \left(T_2\right)={\beta }_0+\frac{{\kappa }_I}{2{\kappa }_R}\ln \left(1+\frac{{\kappa }_R}{2}{T}_2{a_0}^2\right)$

(33)

where a₀ and β₀ are the initial amplitude and phase, respectively. By substituting T₂ = ε²t* and using the Taylor expansion of the phase about ε = 0, the phase is rewritten as equation (34):

$\beta \left(T_2\right)={\beta }_0+\frac{1}{4}{\kappa }_I{T}_0{a_0}^2{\epsilon }^2$

(34)

By substituting equation (34) into the homogeneous solution of equations (i.e., equation (21)), the final format of the homogeneous solution and nonlinear frequency is determined as follows:

$\begin{array}{l}\left\{y_0\right\}=\pm A_{20}{\displaystyle \sum }_{j=1}^8{\left\{A^{\prime }\right\}}_j\hspace{0.5em}{e}^{{\alpha }_{j}X}{e}^{\hspace{0.17em}i{\omega }_{NL}{T}_0}+C.C.\\ A_{20}=\pm \frac{1}{2}\frac{a_0}{\sqrt{1+\frac{{\kappa }_R}{2}{T}_2{a_0}^2}}.e^{i{\beta }_0};\hspace{1em}{\omega }_{NL}=\omega +{\kappa }_{NL}{a_0}^2{\epsilon }^2;\hspace{1em}{\kappa }_{NL}=\frac{1}{4}{\kappa }_I\end{array}$

(35)

where κ _NL is the coefficient of the nonlinear frequency.

4. Results and discussion

To study the accuracy and performance of the presented mathematical formulation, a code has been provided in the Maple environment. In this study, four different types of porosity are discussed. To define the material properties, the following parameters are defined:

$e_0=1-{\alpha }_0^2;\hspace{1em}{e}_m=1-{\alpha }_0$

(36)

where e₀ and e _m are the porosity coefficient and density parameter, respectively. α₀ is the density ratio that is defined as the ratio of the density of the porous material to the density of the dense material. The mechanical properties and volume functions for each type have been reported in Table 1.

Table 1.

Different types of porosity functions.

Porosity type	λ(z)	Г(z)
Power distribution	${\left({\alpha }_0+e_m{\left(\frac{z}{h}+0.5\right)}^g\right)}^2$	${\alpha }_0+e_m{\left(\frac{z}{h}+0.5\right)}^g$
Symmetric distribution	$1-e_0\hspace{0.17em}\cos \left(\frac{\pi z}{h}\right)$	$1-e_m\left(\cos \left(\frac{\pi z}{h}\right)\right)$
Asymmetric distribution	$1-e_0\hspace{0.17em}\cos \left(\frac{\pi z}{2h}+\frac{\pi }{4}\right)$	$1-e_m\left(\cos \left(\frac{\pi z}{2h}+\frac{\pi }{4}\right)\right)$
Uniform distribution	${\left(\frac{2}{\pi }\left(1-{\mathit{e}}_m\right)\right)}^2$	$\sqrt{1-\frac{e_m}{e_0}\left(1-{\left(\frac{2}{\pi }\left(1-e_m\right)\right)}^2\right)}$

where g is the power index for the power distribution porosity, and it is considered one in this parametric study. The inputs for the mathematical code are based on Table 2 and all the reported results are according to these inputs and for simply supported boundary conditions, except those mentioned. Also, the frequencies are in the form of dimensionless and reported values relate to the fundamental frequency.

Table 2.

Material properties for FGP cylinder.

Property (unit)	Quantity	Property (unit)	Quantity
Length (m)	L = 0.8	Maximum bulk modulus (GPa)	K0 = 166.6
Mid-plane radius (m)	R m = 0.16	Maximum shear modulus (GPa)	G0 = 76.9
Thickness (mm)	h = 15	Maximum density (kg/m3)	ρ s = 7800
Poisson’s ratio	ν = 0.3	Density ratio	α0 = 0.65

Table 3 reports the dimensionless linear natural frequencies for homogeneous cylinders versus different L/R _m . In this study, the linear natural frequency means the frequency obtained from the order-zero solution. The results are determined by Eipakchi et al. (2020) (based on classical plate theory), Rao (2007), Eipakchi and Mahboubi Nasrekani (2021) (based on FSDT and used the FE method), and the presented method. To find the results for a homogeneous cylinder using the presented formulation, it is required to assume α₀ = 1. It is seen that the results of the presented method are in good agreement with the other references. The presented method has better accuracy with respect to the FE results in comparison to other references.

Table 3.

Comparison of dimensionless linear natural frequency of homogeneous cylinders (E = 5.096 MPa, R _m /h = 16, R _m = 0.16 m).

L/R m	Present method	FE method (Eipakchi and Mahboubi Nasrekani, 2021)	(Rao, 2007)	(Eipakchi et al., 2020)
2.5	2.635	2.680	2.743	2.616
5	3.370	3.392	3.608	3.435
10	3.431	3.439	3.766	3.572
15	3.438	3.443	3.869	3.592

Table 4 compares the dimensionless linear natural frequencies of FGP cylinders determined by the presented method and the results determined by Daemi and Eipakchi (2019); Daemi and Eipakchi (2019) presented an FE solution for the linear free vibrations of FGP cylinders. The results of the presented method are in good agreement and a maximum error of about 5% is observed with respect to the FE results.

Table 4.

Comparison of dimensionless linear frequency of FGP cylinders (E = 70 GPa, R _m /h = 10, L/R _m = 5, R _m = 0.155 m, ρ _s = 2707).

Porosity type	Presented method	FE (Daemi and Eipakchi, 2019)
Homogeneous (α0 = 1)	3.373	3.204
Symmetric (α0 = 0.89)	3.263	3.056
Asymmetric (α0 = 0.89)	3.255	3.116

Table 5 reports the effect of the boundary conditions on the linear frequencies and the coefficient of nonlinear frequencies. The boundary conditions are shown by two letters, the first letter shows the boundary conditions at x* = 0 and the second letter relates to the boundary conditions at x* = 1 (e.g., S-C means simply support at x* = 0 and clamped at x* = 1). The nonlinear frequency behavior is shown by two letters where H stands for hardening behavior and S stands for softening behavior. It is seen that for all boundary conditions except C-F, the nonlinear frequency behavior is hardening, and the nonlinear frequency increases in comparison to the linear frequency for the studied cases. It is also concluded that by using porous materials we can decrease the weight of the cylinder by about 22% (for the symmetric distribution) while the frequency drops by about 9% with respect to the homogeneous cylinder in the absence of porosity. For all the cases of porosity, the maximum and minimum frequencies relate to the C-C and S-F boundary conditions, respectively. These findings have significant implications for the design and optimization of lightweight, vibration-resistant structures. By strategically incorporating porosity and tailoring boundary conditions, engineers can fine-tune the dynamic response of cylindrical shells to meet specific performance requirements.

Table 5.

Boundary condition effects on FGP cylinder linear and nonlinear frequencies.

Porosity type		Homogeneous (α0 = 1)	Power distribution	Symmetric	Asymmetric	Uniform
Mass (kg)		94.0	77.8	73.1	73.4	66.4
C-C	ω	3.401	3.114	3.070	3.070	1.675
	κ NL	1.332	1.080	1.805	1.067	0.656
	Behavior (H/S)	H	H	H	H	H
S-S	ω	3.373	3.056	3.043	3.007	1.661
	κ NL	2.125	3.320	2.933	3.991	1.005
	Behavior (H/S)	H	H	H	H	H
C-S	ω	3.387	3.085	3.056	3.038	1.668
	κ NL	1.532	0.586	2.345	1.063	0.754
	Behavior (H/S)	H	H	H	H	H
C-F	ω	1.721	1.577	1.553	1.555	0.847
	κ NL	−0.022	−0.022	−0.014	−0.022	−0.011
	Behavior (H/S)	S	S	S	S	S
S-F	ω	1.716	1.565	1.547	1.541	0.845
	κ NL	0.978	0.082	0.673	0.091	0.482
	Behavior (H/S)	H	H	H	H	H

Figure 2 shows the effect of the density ratio (α₀) on the mass and linear frequency for different types of FGP cylinders. It is seen that in the studied range of parameters, for the symmetric FGP cylinder, we can decrease the weight of the structure by about 60% by decreasing the density ratio while the linear frequency remains nearly the same.

Figure 2.

Effects of density ratio on total mass and dimensionless linear frequency of different FGP cylinder.

Figures 3 and 4 show the effect of density ratio on the nonlinear behavior of different FGP cylinders with different types of porosity. It is seen that only for the power and asymmetric distributions when α₀ = 0.1, the nonlinear frequency behavior is softening and for other cases, it is hardening. For all types of porosity distribution, the maximum variation of nonlinear frequency occurs when the density ratio is about 0.1 for the studied range. It can be concluded that for the studied range of parameters, the density ratio of about 0.3 is an optimum value. Because, according to the results of Figure 2, the mass is comparatively low, and nonlinear frequency variation is negligible (except for the uniform distribution).

Figure 3.

Effect of density ratio on nonlinear frequency for different FGP cylinders with symmetric and asymmetric distributions.

Figure 4.

Effect of density ratio on nonlinear frequency for different FGP cylinders with uniform and power distributions.

Figure 5 shows the effect of R _m /h ratio on the nonlinear frequency for an FGP cylinder with symmetric porosity distribution. It is seen that by increasing R _m /h the nonlinear frequency behavior of the FGP cylinder converts from hardening to softening. It should be noted that Rm/h ratio is a parameter to categorize the cylinder into thin or thick groups. It means that by increasing Rm/h ratio, the cylinder can be categorized into thin shells and the stress distribution across the thickness can be assumed uniformly. It might be the reason for the change in the nonlinear frequencies from hardening to softening.

Figure 5.

Effect of R _m /h ratio on nonlinear frequency of FGP cylinder with symmetric porosity distribution.

Figure 6 reports the effect of the L/R _m ratio on the nonlinear frequency for an FGP cylinder with symmetric porosity distribution. For all studied cases, the nonlinear frequency behavior is hardening. The same results are observed for other types of porosity distribution which are not reported here.

Figure 6.

Effect of L/R _m ratio on nonlinear frequency of FGP cylinder with symmetric porosity distribution.

5. Conclusion

In this study, an analytical solution based on the multiple-scale method was presented to analyze the nonlinear free vibration behavior of FGP cylinders. To extract the governing equations, Mirsky–Hermann’s theory was employed and to consider the moderately large deformations, the von Karman relations were used. By performing a parametric study, the effect of different parameters such as porosity distributions, boundary conditions, and geometrical parameters were studied on the frequencies. In the studied range of materials and geometric dimensions:

• For all boundary conditions except C-F, the nonlinear frequency behavior is hardening.

• Existence of porosity decreases the weight of the cylinder by about 60% with respect to non-porous cylinder while the frequency remains nearly the same.

• The maximum and minimum frequencies occur at the C-C and S-F boundary conditions, respectively.

• Only for the power and asymmetric distributions when α₀ = 0.1, the nonlinear frequency behavior is softening and for other cases, it is hardening.

• For all types of studied porosity distributions, when the density ratio is about 0.1, the variation of the nonlinear frequency versus amplitude is significant.

• By increasing R _m /h ratio the nonlinear frequency behavior converts from hardening to softening.