Analysis of Geometrically Nonlinear Vibrations of Functionally Graded Shallow Shells of a Complex Shape

Abstract

Geometrically nonlinear vibrations of functionally graded shallow shells of complex planform are studied. The paper deals with a power-law distribution of the volume fraction of ceramics and metal through the thickness. The analysis is performed with the use of the R-functions theory and variational Ritz method. Moreover, the Bubnov-Galerkin and the Runge-Kutta methods are employed. A novel approach of discretization of the equation of motion with respect to time is proposed. According to the developed approach, the eigenfunctions of the linear vibration problem and some auxiliary functions are appropriately matched to fit unknown functions of the input nonlinear problem. Application of the R-functions theory on every step has allowed the extension of the proposed approach to study shallow shells with an arbitrary shape and different kinds of boundary conditions. Numerical realization of the proposed method is performed only for one-mode approximation with respect to time. Simultaneously, the developed method is validated by investigating test problems for shallow shells with rectangular and elliptical planforms, and then applied to new kinds of dynamic problems for shallow shells having complex planforms.

Keywords:
Functionally graded shallow shells; R-functions theory; numerical-analytical approach; complex planform

1 INTRODUCTION

Structural elements modeled by shallow shells are widely used in various engineering fields, including, for instance: mechanical, aerospace, marine, military, and civil engineering. Such elements can have various planforms, boundary conditions (including mixed ones), and different types of curvature. In order to improve the strength of the modern design, a new class of composite materials, i.e., functionally graded materials (FGM), has been recently applied. In spite of the observation that FGM are inhomogeneous, they have essential advantage, i.e., material properties undergo smooth and continuous variations in the thickness direction. This is why the stress concentration present in laminated structures can be eliminated. However, an analysis of functionally graded (FG) shells is more complicated than of homogeneous material structures, since partial differential equations (PDEs) with variable coefficients govern the shallow shells made of FGM, and the strain-stress fields are coupled. As it is known, getting a validated solution to the mentioned PDEs creates a very difficult problem even in the case of shells having relatively simple planforms. In addition, the problem becomes more complicated if FG shallow shells perform vibrations at large amplitudes.

This class of problems is challenging and there exist numerous investigations devoted to analysis of dynamical behavior of the FG plates and shells. This is especially true for linear problems (Lam et al. (2002); Loy et al. (2008); Matsunaga (2008); Neves et al. (2013); Pradyumna and Bandyopadhyay (2008); Reddy et al. (1999); Reddy (2009); Shen (2009); Tornabene and Viola (2007, 2009); Tornabene et al. (2011); Zhu et al. (2014, 2015)). However, in the last decade, nonlinear free and forced vibrations of the FG shells have been extensively studied as well (Reddy et al. (1999); Reddy (2009); Shen (2009); Amabili (2008); Alijani et al. (2011a,b); Bich et al. (2012); Chorfi and Houmat (2010); Sundararajan et al. (2005); Woo et al. (2006); Xiang et al. (2015); Zhao and Liew (2009)). A complete survey on linear and nonlinear vibrations of FG plates and shells can be found in the following references (Shen (2009); Tornabene et al. (2011); Zhu et al. (2014); Amabili (2008); Alijani et al. (2011b)). Note that in the aforementioned papers, nonlinear vibrations of simply supported or clamped FG structures of rectangular, skew or circle planforms have been analyzed by different numerical methods such as Finite Element Method (FEM) [Lam et al. (2002); Loy et al. (2008); Matsunaga (2008); Reddy et al. (1999); Reddy (2009); Shen (2009); Bich et al. (2012)], Differential Quadrature Method [Tornabene and Viola (2007, 2009); Tornabene et al. (2011)], modified Fourier-Ritz Approach [Zhu et al. (2014, 2015)], etc. A survey on vibrations of open shells of revolution can be found in papers [Zhu et al. (2014, 2015); Amabili (2008)]. The analysis of published literature devoted to nonlinear vibrations of FG shallow shells is, in general, restricted to their simple planforms and classical boundary conditions. However, FG shells of arbitrary planforms and different boundary conditions are widely used in practice. Consequently, it is important to develop universal and effective methods for investigation of nonlinear vibrations of functionally graded shallow shells of complex planforms and different boundary conditions. Earlier, in papers [Kurpa (2009a); Kurpa et al. (2007, 2010, 2013); Awrejcewicz et al. (2010, 2013, 2015a); Kurpa and Mazur (2010); Kurpa and Shmatko (2014)], the original meshless method aimed at an application of the R-functions theory, variational Ritz method, Bubnov-Galerkin procedure, and Runge-Kutta method has been proposed. In reference [Awrejcewicz et al. (2015b)] this method has been extended to geometrically nonlinear vibration problems of functionally graded shallow shells of arbitrary planforms. In the aforementioned paper, theoretical results have been presented in a rather short form.

In this paper we present a numerical-analytical method based on the R-functions theory in more detail. Formulation of the problem is performed using theories of the shallow shells, i.e., classical (CTS) and geometrically refined nonlinear theory of the first order (FSDT). Distinctive feature of the proposed approach is the original method of reducing the input nonlinear system of PDEs to a nonlinear system of ordinary differential equations. In order to realize this approach appropriately, it is needed to perform a free vibration analysis of FG shells with high accuracy, since eigenfunctions of the linear vibration problem are used on each step of the proposed algorithm. The proposed method is validated by investigating test problems for shallow shells of rectangular and elliptical planforms, and then applied to new vibration problems for shallow shells of complex planforms.

2 MATHEMATICAL FORMULATION

Consider a functionally graded shallow shell with uniform thickness h, made of a mixture of ceramics and metal. It is assumed that the shell has an arbitrary planform. The volume fraction of ceramic _^Vc and metal phases _^Vm are related by formula _{^{Vc + Vm =}} 1. The volume fraction of ceramics _^Vc can be expressed as

$V_c={\left(\frac{z}{h}+\frac{1}{2}\right)}^k.$ (1)

In the above formula, the index k (0 ≤ k < ∞) denotes the volume fraction exponent, z is the distance between the current point and the shell midsurface. In the case when the power index k is equal to zero, one obtains a homogeneous material (ceramics), but if k approaches infinity, the shell is purely metallic.

It should be noted that FG materials are widely used in high-temperature environments and their mechanical characteristics can be different, depending on temperature changes. Therefore, this temperature dependence must be taken into account to obtain more accurate solution. We use the relations given in [Reddy (2009); Shen (2009)] in the following form

$P_j(T)=P_0\left(P_{-1}{T}^{-1}+1+P_{1}T+P_2{T}^2+P_3{T}^3\right),$

where P ₀, P _-1 P ₁ P ₂ P ₃ are coefficients determined for each specific material. A table of values of these coefficients for some materials is presented in [Reddy (2009); Shen (2009); Alijani et al. (2011a)].

As it is known [Reddy et al. (1999); Reddy (2009); Shen (2009); Tornabene and Viola (2007, 2009); Tornabene et al. (2011); Alijani et al. (2011a,b); Bich et al. (2012)], in FGM structures the material properties along the thickness are proportional to the volume fraction of the constituent materials, i.e., we have

$P\left(z,T\right)=\left(P_c\left(T\right)-P_m\left(T\right)\right){\left(\frac{z}{h}+\frac{1}{2}\right)}^k+P_m\left(T\right),$

or equivalently,

$P\left(z,T\right)=\left(P_c\left(T\right)-P_m\left(T\right)\right)V_c+P_m\left(T\right),$ (2)

where _^Pc (T), _^Pm (T) are the corresponding characteristics of the ceramic and metal, respectively. Equation (2) represents a general formula for determination of the elastic modulus E, Poisson's ratio v, and the density ρ of the composite, that is

$E=\left(E_c-E_m\right)V_c+E_m,\nu =\left({\nu }_c-{\nu }_m\right)V_c+\nu E_m,\rho =\left({\rho }_c-{\rho }_m\right)V_c+{\rho }_m.$ (3)

According to the nonlinear first-order shear deformation theory of shallow shells (FSDT), displacement components _^ui , u ₂, u ₃ at a point (x, y, z) are expressed as functions of the shell middle surface displacements u, v, and w in the Ox, Oy, and Oz directions, and the independent rotations _{^{ψx, ψy}} of the transverse normal to the middle surface about the Oy and Ox axes, respectively as [Reddy (2009); Reddy et al. (1999)]:

$u_1=u+z{\psi }_x,u_2=v+z{\psi }_y,u_3=w.$

Strains ε = {ε₁₁; ε₂₂; ε₁₂}^T at an arbitrary point of the shallow shell follow

${\epsilon }_{11}={\epsilon }_{11}^0+z{\chi }_{11},{\epsilon }_{22}={\epsilon }_{22}^0+z{\chi }_{22},{\epsilon }_{12}={\epsilon }_{12}^0+z{\chi }_{12},$ (4)

where

${\epsilon }_{11}^0=\frac{\partial u}{\partial x}-k_{1}w+\frac{1}{2}{\left(\frac{\partial w}{\partial x}\right)}^2,{\epsilon }_{22}^0=\frac{\partial v}{\partial y}-k_{2}w+\frac{1}{2}{\left(\frac{\partial w}{\partial y}\right)}^2,{\epsilon }_{12}^0=\frac{\partial u}{\partial y}+\frac{\partial v}{\partial x}+\frac{\partial w}{\partial x}\frac{\partial w}{\partial y},$ (5)

and k ₁ = 1 / ^_Rx , k ₂ = 1/^_Ry are the principal curvatures of the shell along the coordinates x and y, respectively.

Let us present formulas (5) employing the following general notation:

$\left\{\epsilon \right\}=\left\{{\epsilon }^{L_0}\right\}+\left\{{\epsilon }^N\right\},$ (6)

where

$\left\{{\epsilon }^{L_0}\right\}={\left\{{\epsilon }_{11}^{L_0},{\epsilon }_{22}^{L_0},{\epsilon }_{12}^{L_0}\right\}}^T={\left\{\frac{\partial u}{\partial x}-k_{1}w;\frac{\partial v}{\partial y}-k_{2}w;\left(\frac{\partial u}{\partial y}+\frac{\partial v}{\partial x}\right)\right\}}^T,$ (7)

$\begin{array}{l}\left\{{\epsilon }^{N)}\right\}={\left\{{\epsilon }_{11}^N,{\epsilon }_{22}^N,{\epsilon }_{12}^N\right\}}^T={\left\{\frac{1}{2}{\left(\frac{\partial w}{\partial x}\right)}^2;\frac{1}{2}{\left(\frac{\partial w}{\partial y}\right)}^2;\frac{\partial w}{\partial x}\cdot \frac{\partial w}{\partial y}\right\}}^T,\\ {\epsilon }_{13}=\delta \left(\frac{\partial w}{\partial x}+{\psi }_x\right),{\epsilon }_{23}=\delta \left(\frac{\partial w}{\partial y}+{\psi }_y\right).\end{array}$ (8)

Strains χ = {χ ₁₁; χ ₂₂; χ ₁₂}^T are defined by the following formulas

${\chi }_{11}=\delta \frac{\partial {\psi }_x}{\partial x}-\left(1-\delta \right)\frac{{\partial }^{2}w}{\partial x^2},{\chi }_{22}=\delta \frac{\partial {\psi }_y}{\partial y}-\left(1-\delta \right)\frac{{\partial }^{2}w}{\partial y^2},{\chi }_{12}=\delta \left(\frac{\partial {\psi }_x}{\partial y}+\frac{\partial {\psi }_y}{\partial x}\right)-\left(1-\delta \right)\frac{{\partial }^{2}w}{\partial x\partial y},$ (9)

where indicator δ is the tracing constant taking values 1 or 0 for the FSDT and CST, respectively.

The strain resultants N = (N ₁₁, N ₂₂, N ₁₂)^T moment resultants M = (M ₁₁, M ₂₂, M ₁₂)^T and shear stress resultants Q = (^_Qx , ^_Qy )^T are calculated by integration along Oz-axis, and they have the following forms

$N=\left[A\right]\left\{\epsilon \right\}+\left[B\right]\left\{\chi \right\},M=\left[B\right]\left\{\epsilon \right\}+\left[D\right]\left\{\chi \right\},$ (10)

$\left[A\right]=\left[\begin{array}{l}{A}_{11}{A}_{12}0\\ A_{12}{A}_{22}0\\ 00A_{33}\end{array}\right],\left[B\right]=\left[\begin{array}{l}{B}_{11}{B}_{12}0\\ B_{12}{B}_{22}0\\ 00B_{33}\end{array}\right],\left[D\right]=\left[\begin{array}{l}{D}_{11}{D}_{12}0\\ D_{12}{D}_{22}0\\ 00D_{33}\end{array}\right],$ (11)

where

$\left(\left[A\right],\left[B\right],\left[D\right]\right)={\displaystyle \underset{-\frac{h}{2}}{\overset{\frac{h}{2}}{\int }}Q\left(z\right)\left(\mathrm{1,}z,z^2\right)dz},Q\left(z\right)=\frac{E\left(z\right)}{1-{\nu }^2\left(z\right)}\left[C\right],$ (12)

$\left[C\right]=\left[\begin{array}{l}1\nu 0\\ \nu 10\\ 00\frac{1-\nu }{2}\end{array}\right].$ (13)

Transverse shear force resultants _^Qx , _^Qy are defined as

$Q_x=K_s^2{A}_{33}{\epsilon }_{13},Q_y=K_s^2{A}_{33}{\epsilon }_{23},$ (14)

where $K_s^2$ denotes the shear correction factor (in this paper it is selected to be 5/6).

Further, let us consider materials with Poisson’s ratio independent of the temperature, being the same for ceramic and metal, i.e., _^vm = _^vc . In this case, the in-plane force resultants N = (N ₁₁, N ₂₂, N ₁₂)^T and moment resultants M = (M ₁₁, M ₂₂, M ₁₂)^T in the framework of FSDT, taking into account the power law (1), are defined as follows

$\left\{N\right\}=\frac{1}{1-{\nu }^2}\left[C\right]\left(E_1\left\{\epsilon \right\}+E_2\left\{\chi \right\}\right),\left\{M\right\}=\frac{1}{1-{\nu }^2}\left[C\right]\left(E_1\left\{\epsilon \right\}+E_3\left\{\chi \right\}\right),$ (15)

$E_1=\left(E_m+\frac{E_c-E_m}{k+1}\right)h,E_2=\frac{\left(E_c-E_m\right)k{h}^2}{2(k+1)\left(k+2\right)},$ (16)

$E_3=\left(\frac{E_m}{12}+\left(E_c-E_m\right)\left(\frac{1}{k+3}-\frac{1}{k+2}+\frac{1}{4(k+1)}\right)\right)h^3.$ (17)

Mass density ρ is also estimated by integration along the shell thickness, to yield

$\rho =\left({\rho }_m+\frac{{\rho }_c-{\rho }_m}{k+1}\right)h.$ (18)

3 SOLUTION METHOD

3.1 Linear Problem

The first step of the proposed method is to solve the linear vibration problem. For this purpose, the vector of unknown functions is represented as

$\begin{array}{l}\overrightarrow{U}(\overline{u}(x,y,t),\overline{v}(x,y,t),\overline{w}(x,y,t),{\overline{\psi }}_x(x,y,t),{\overline{\psi }}_y(x,y,t))=\\ =\overrightarrow{U}(u(x,y,t),v(x,y,t),w(x,y,t),{\psi }_x(x,y,t),{\psi }_y(x,y,t)),\end{array}$

where λ is a vibration frequency. Applying Hamilton’s principle, we get a variational equation in the following form

$\partial \left(U_{\max }-\lambda T_{\max }\right)=\mathrm{0,}$ (19)

where U and T are strain and potential energy, respectively, and

$\begin{array}{l}{U}_{\max }=\frac{1}{2}\underset{\Omega }{\int \int }\left(N_{11}^L{\epsilon }_{11}^{L_0}+N_{22}^L{\epsilon }_{22}^{L_0}+N_{12}^L{\epsilon }_{12}^{L_0}+M_{11\chi 11}^L+M_{22\chi 22}^L+M_{12\chi 12}^L+\delta \left(Q_x{\epsilon }_{13}+Q_y{\epsilon }_{23}\right)\right)d\Omega ,\\ T_{\max }=\frac{1}{2}\underset{\Omega }{\int \int }\left(I_0\left(u^2+v^2+w^2\right)+2I_1\delta \left(u{\psi }_x+v{\psi }_y\right)+I_2\delta \left({\psi }_x^2+{\psi }_y^2\right)\right)d\Omega ,\end{array}$

whereas $\left\{{\epsilon }^{L_0}\right\}=\left\{{\epsilon }_{11}^{L_0},{\epsilon }_{22}^{L_0},{\epsilon }_{12}^{L_0}\right\}$ and {χ} = {χ₁₁, χ₂₂, χ₁₂} are defined by formulas (7), (9), and

$\begin{array}{l}\left\{N^L\right\}={\left\{N_{11}^L,N_{22}^L,N_{12}^L\right\}}^T=\frac{1}{1-{\nu }^2}\left[C\right]\left(E_1\left\{{\epsilon }^{L_0}\right\}+E_2\left\{\chi \right\}\right),\\ \left\{M^L\right\}={\left\{M_{11}^L,M_{22}^L,M_{12}^L\right\}}^T=\frac{1}{1-{\nu }^2}\left[C\right]\left(E_1\left\{\epsilon \right\}+E_3\left\{\chi \right\}\right),\\ I_0=\left({\rho }_m+\frac{{\rho }_c-{\rho }_m}{k+1}\right)h,\begin{array}{cc}& \end{array}{I}_1{=}_{-\frac{h}{2}}^{\int \frac{h}{2}}\rho \left(z\right)zdz=\frac{k\left({\rho }_c-{\rho }_m\right)}{2\left(k+1\right)\left(k+2\right)}{h}^2,\\ I_2{=}_{-\frac{h}{2}}^{\int \frac{h}{2}}\rho \left(z\right)z^{2}dz=\left(\frac{{\rho }_m}{12}+\left({\rho }_c-{\rho }_m\right)\left(\frac{1}{k+3}-\frac{1}{k+2}+\frac{1}{4\left(k+1\right)}\right)\right)h^3.\end{array}$ (20)

Minimization of functional (19) is performed using the Ritz method, and the necessary sequence of coordinate functions is built with the help of the R-functions theory [Rvachev (1982); Kurpa (2009b)].

3.2 Nonlinear Problem

For simplicity, let us describe the algorithm for solution to the nonlinear vibration problem in the frame of the classical theory (CTS). Note that inertia terms are ignored in motion equations while solving the nonlinear problem. Let us formulate the given problem with respect to shell displacements:

$L_{11}u+L_{12}v+L_{13}w=N{L}_1\left(w\right),$ (21)

$L_{21}u+L_{22}v+L_{23}w=N{L}_2\left(w\right),$ (22)

$L_{31}u+L_{32}v+L_{33}w=N{L}_{32}\left(u,v,w\right)+N{L}_{33}\left(w\right)+m_1\frac{{\partial }^{2}w}{\partial t^2},$ (23)

where linear differential operators _^Lij (i, j, = 1, 2, 3) are defined in the following way

$\begin{array}{l}{L}_{11}=A_{11}\frac{{\partial }^2}{\partial x^2}+A_{33}\frac{{\partial }^2}{\partial y^2},L_{12}=L_{21}=\left(A_{12}+A_{33}\right)\frac{{\partial }^2}{\partial x\partial y},\\ L_{13}=-L_{31}=-B_{11}\frac{{\partial }^3}{\partial x^3}-\left(B_{12}+2B_{33}\right)\frac{{\partial }^3}{\partial x\partial y^2}-\left(k_1{A}_{11}+k_2{A}_{12}\right)\frac{\partial }{\partial x},\\ L_{22}=A_{22}\frac{{\partial }^2}{\partial y^2}+A_{33}\frac{{\partial }^2}{\partial x^2},L_{23}-L_{32}=B_{22}\frac{{\partial }^3}{\partial y^3}-\left(B_{12}+2B_{33}\right)\frac{{\partial }^3}{\partial x^2\partial y}-\left(k_1{A}_{21}+k_2{A}_{22}\right)\frac{\partial }{\partial y},\\ L_{13}=L_{13}^{\left(pl\right)}+L_{13}^{\left(cur\right)},L_{23}=L_{23}^{\left(pl\right)}+L_{23}^{\left(cur\right)},L_{33}=L_{33}^{\left(pl\right)}+L_{33}^{\left(cur\right)},\\ L_{13}^{\left(pl\right)}=-B_{11}\frac{{\partial }^3}{\partial x^3}-\left(B_{12}+2B_{33}\right)\frac{{\partial }^3}{\partial x\partial y^2},L_{13}^{\left(cur\right)}=-\left(k_1{A}_{11}+k_2{A}_{12}\right)\frac{\partial }{\partial x},\\ L_{23}^{\left(pl\right)}=-B_{22}\frac{{\partial }^3}{\partial y^3}-\left(B_{21}+2B_{33}\right)\frac{{\partial }^3}{\partial x^2\partial y},L_{23}^{\left(cur\right)}=-\left(k_1{A}_{21}+k_{2}A22\right)\frac{\partial }{\partial y},\\ L_{33}^{\left(pl\right)}=-\left(D_{11}\frac{{\partial }^4}{\partial x^4}+2\left(D_{12}+2D_{33}\right)\frac{{\partial }^4}{\partial x^2\partial y^2}+D_{22}\frac{{\partial }^4}{\partial y^4}\right),\\ L_{33}^{\left(cur\right)}=-2\left(k_1{B}_{11}+k_2{B}_{12}\right)\frac{{\partial }^2}{\partial x^2}-2\left(k_1{B}_{21}+k_2{B}_{22}\right)\frac{{\partial }^2}{\partial y^2}-\left(k_1^2{A}_{11}+2k_1{k}_2{A}_{12}+k_2^2{A}_{22}\right).\end{array}$

The expressions of the nonlinear terms standing on the right-hand side of the system (21)-(23) are defined as follows

$\begin{array}{l}N{L}_1\left(w\right)=-L_{11}\left(w\right)\frac{\partial w}{\partial x}-L_{12}\left(w\right)\frac{\partial w}{\partial y},N{L}_2\left(w\right)=-L_{12}\left(w\right)\frac{\partial w}{\partial x}-L_{22}\left(w\right)\frac{\partial w}{\partial y},\\ N{L}_{32}\left(u,v,w\right)=-N_{11}^L\frac{{\partial }^{2}w}{\partial x^2}+2N_{12}^L\frac{{\partial }^{2}w}{\partial x\partial y}+N_{22}^L\frac{{\partial }^{2}w}{\partial y^2}-\left(\frac{{\partial }^2{M}_{11}^N}{\partial x^2}+\frac{{\partial }^2{M}_{22}^N}{\partial y^2}+\frac{{\partial }^2{M}_{12}^N}{\partial x\partial y}\right)-\\ -\frac{\partial w}{\partial x}\left(L_{31}^{\left(pl\right)}+\frac{1}{2}{L}_{31}^{\left(cur\right)}\right)w-\frac{\partial w}{\partial y}\left(L_{32}^{\left(pl\right)}+\frac{1}{2}{L}_{32}^{\left(cur\right)}\right)w,\\ N{L}_{33}\left(w\right)=-N_{11}^{NL}\frac{{\partial }^{2}w}{\partial x^2}-2N_{12}^{NL}\frac{{\partial }^{2}w}{\partial x\partial y}-N_{22}^{NL}\frac{{\partial }^{2}w}{\partial y^2},\end{array}$

where: $\left\{N^L\right\}={\left\{N_{11}^L,N_{22}^L,N_{12}^L\right\}}^T,\left\{M^L\right\}={\left\{M_{11}^L,M_{22}^L,M_{12}^L\right\}}^T$ are defined by formulas (20), and $\left\{N^{NL}\right\}={\left\{N_{11}^{NL};N_{22}^{NL};N_{12}^{NL}\right\}}^T=\left[A\right]\left\{{\epsilon }^N\right\},\left\{M^N\right\}={\left\{M_{11}^N;M_{22}^N;M_{12}^{NL}\right\}}^T=\left[B\right]\left\{{\epsilon }^N\right\}$.

The vector {_^(N } is defined by formula (8).

Note that the type of linear operators L ₁₃, L ₂₃, L ₃₃ is simplified in the case of plates due to the condition k ₁ = k ₂ = 0, which implies $L_{13}^{\left(cur\right)}=L_{23}^{\left(cur\right)}=L_{33}^{\left(cur\right)}$ = 0.

Let us take the unknown functions w(x, y, t) u(x, y, t), v(x, y, t) in the form of an expansion in terms of eigenfunctions $w_i^{\left(e\right)}\left(x,y\right),u_i^{\left(e\right)}\left(x,y\right),v_i^{\left(e\right)}\left(x,y\right)$ of the linear problem with coefficients _^yi (t) depending on time, i.e. we have

$w={\displaystyle \sum _{i=1}^n{y}_i}\left(t\right)w_i^{\left(e\right)}\left(x,y\right),$ (24)

$u={\displaystyle \sum _{i=1}^n{y}_i}\left(t\right)u_i^{\left(e\right)}\left(x,y\right)+{\displaystyle \sum _{i=1}^n{\displaystyle \sum _{j=1}^n{y}_i{y}_j{u}_{ij},\begin{array}{c}\end{array}}}v={\displaystyle \sum _{i=1}^n{y}_i\left(t\right)v_i^{\left(e\right)}\left(x,y\right)+{\displaystyle \sum _{i=1}^n{\displaystyle \sum _{j=1}^n{y}_i{y}_j{v}_{ij}.}}}$ (25)

The functions _^uij , _^vij should satisfy the following system of differential equations

$\{\begin{array}{l}{L}_{11}\left(u_{ij}\right)+L_{12}\left(v_{ij}\right)=-N{L}_1^{\left(2\right)}\left(w_i^{\left(e\right)},w_j^{\left(e\right)}\right),\\ L_{21}\left(u_{ij}\right)+L_{22}\left(v_{ij}\right)=-N{L}_2^{\left(2\right)}\left(w_i^{\left(e\right)},w_j^{\left(e\right)}\right)\end{array}$ (26)

where

$\begin{array}{l}N{L}_1^{(2)}\left(w_i^{(e)},w_j^{(e)}\right)=w_i^{\left(e\right)}{,}_x{L}_{11}{w}_j^{\left(e\right)}+w_i^{\left(e\right)}{,}_y{L}_{12}{w}_j^{\left(e\right)},\\ N{L}_2^{(2)}\left(w_i^{(e)},w_j^{(e)}\right)=w_i^{\left(e\right)}{,}_x{L}_{12}{w}_j^{\left(e\right)}+w_i^{\left(e\right)}{,}_y{L}_{22}{w}_j^{\left(e\right)}.\end{array}$

The system of equations (26) can be solved with the RFM for any planform and various kinds of boundary conditions. It is possible to show that the variational formulation of this problem is reduced to finding the minimum of the following functional:

$\begin{array}{l}I\left(u_{ij},v_{ij}\right)=\underset{\Omega }{{\displaystyle \int \int }}\lfloor \left(N_{11}^{\left(L_2\right)}{\epsilon }_{11}^{L_2}+N_{22}^{\left(L_2\right)}{\epsilon }_{22}^{L_2}+N_{12}^{\left(L_2\right)}{\epsilon }_{12}^{L_2}\right)-\\ -2\left(N{L}_1\left(w_i^{\left(e\right)},w_j^{\left(e\right)}\right)u_{ij}+N{L}_2\left(w_i^{\left(e\right)},w_j^{\left(e\right)}\right)v_{ij}\right)\rfloor d\Omega -2\underset{\partial \Omega }{\int }\left(F_1{u}_{ij}^{\left(n\right)}+F_2{v}_{ij}^{\left(n\right)}\right)ds,\end{array}$ (27)

where $N_{ij}^{\left(L_2\right)}$ (i, j = 1, 2) are defined by the following expressions

$\left\{N^{\left(L_2\right)}\right\}={\left\{N_{11}^{\left(L_2\right)};N_{22}^{\left(L_2\right)};N_{12}^{\left(L_2\right)}\right\}}^T=\left[A\right]{\left\{{\epsilon }^{L_2}\right\}}^T,\begin{array}{cc}& \end{array}\left\{{\epsilon }^{L_2}\right\}=\left\{\frac{\partial u_{ij}}{\partial x};\frac{\partial v_{ij}}{\partial y};\frac{\partial u_{ij}}{\partial y}+\frac{\partial v_{ij}}{\partial x}\right\}.$ (28)

The functions F ₁ and F ₂ are:

$\begin{array}{l}{F}_1=-\left(N_{11}^{\left(N_2\right)}{l}^2+N_{22}^{\left(N_2\right)}{m}^2+2N_{12}^{\left(N_2\right)}lm\right),\\ F_2=-N_{12}^{\left(N_2\right)}\left(l^2-m^2\right)-lm\left(N_{22}^{\left(N_2\right)}-N_{11}^{\left(N_2\right)}\right),\\ \left\{N^{\left(N_2\right)}\right\}={\left\{N_{11}^{\left(N_2\right)};N_{22}^{\left(N_2\right)};N_{12}^{\left(N_2\right)}\right\}}^T=\left[A\right]\left\{{\epsilon }^{N_2}\left(w_i^{\left(e\right)},w_j^{\left(e\right)}\right)\right\},\end{array}$ (29)

$\left\{{\epsilon }^{N_2}\left(w_i^{\left(e\right)},w_j^{\left(e\right)}\right)\right\}=\left\{{\epsilon }_{11}^{N_2},{\epsilon }_{22}^{N_2},{\epsilon }_{12}^{N_2}\right\}=\frac{1}{2}{\left\{\left(\frac{\partial w_i^{\left(e\right)}}{\partial x}\right)\left(\frac{\partial w_j^{\left(e\right)}}{\partial x}\right);\left(\frac{\partial w_i^{\left(e\right)}}{\partial y}\right)\left(\frac{\partial w_j^{\left(e\right)}}{\partial y}\right);\left(\frac{\partial w_i^{\left(e\right)}}{\partial x}\frac{\partial w_j^{\left(e\right)}}{\partial y}+\frac{\partial w_j^{\left(e\right)}}{\partial x}\frac{\partial w_i^{\left(e\right)}}{\partial y}\right)\right\}}^T,$ (30)

where l and m are directional cosines of the normal n to the border.

Observe that in case of boundary conditions with clamped edge, we have

$u_{ij}^{\left(n\right)}=\mathrm{0,}{v}_{ij}^{\left(n\right)}=\mathrm{0,}$

and consequently, a contour integral in formula (27) equals zero. The system of basic functions for functional (27) is built with the help of the R-functions theory.

Substituting expressions (24)-(25) for the functions u, v, w into equations of motion (21)-(23) and applying the Bubnov-Galerkin procedure, we obtain the following system of nonlinear ordinary differential equations for the unknown functions _^yr (t), (r = 1, ..., n):

$y_r^n\left(t\right)+{\omega }_{Lr}^2{y}_r\left(t\right)+{\displaystyle \sum _{i,j=1}^n{\beta }_{ij}^{\left(r\right)}{y}_i\left(t\right)y_j\left(t\right)+{\displaystyle \sum _{i,j,k=1}^n{\gamma }_{ijk}^{\left(r\right)}{y}_i\left(t\right)y_j\left(t\right)y_k\left(t\right)=0.}}$ (31)

The formulas defining the coefficients ${\beta }_{ij}^{\left(r\right)},{\gamma }_{ijk}^{\left(r\right)}$ follow:

$\begin{array}{l}{\beta }_{ij}^{\left(r\right)}=\frac{-1}{m_1{\Vert w_r^{\left(e\right)}\Vert }^2}\underset{\Omega }{\int \int }(N_{11}^L\left(u_i^{\left(e\right)},v_i^{\left(e\right)},w_i^{\left(e\right)}\right)\frac{{\partial }^2{w}_j^{\left(e\right)}}{\partial x^2}+N_{22}^L\left(u_i^{\left(e\right)},v_i^{\left(e\right)},w_i^{\left(e\right)}\right)\frac{{\partial }^2{w}_j^{\left(e\right)}}{\partial y^2}+2N_{12}^L\left(u_i^{\left(e\right)},v_i^{\left(e\right)},w_i^{\left(e\right)}\right)\frac{{\partial }^2{w}_j^{\left(e\right)}}{\partial x\partial y}+\\ +\frac{{\partial }^2{M}_{11}^{N_{2p}}}{\partial x^{}}+\frac{{\partial }^2{M}_{22}^{N_{2p}}}{\partial y^2}+\frac{{\partial }^2{M}_{12}^{N_{2p}}}{\partial x\partial y}-k_1{N}_{11}^{N_{2p}}\frac{\partial w_i^{\left(e\right)}}{\partial x}-k_2{N}_{22}^{N_{2p}}\frac{\partial w_i^{\left(e\right)}}{\partial x})w_r^{\left(e\right)}d\Omega ,\end{array}$ (32)

$\begin{array}{l}{\gamma }_{ijk}^{\left(r\right)}=-\frac{1}{m_1{\Vert w_r^{\left(e\right)}\Vert }^2}\underset{\Omega }{\int \int }(N_{11}^{N_{2p}}\left(u_{ij},v_{ij},w_i^{\left(e\right)},w_j^{\left(e\right)}\right)\frac{{\partial }^2{w}_k^{\left(e\right)}}{\partial x^2}+N_{22}^{N_{2p}}\left(u_{ij},v_{ij},w_i^{\left(e\right)},w_j^{\left(e\right)}\right)\frac{{\partial }^2{w}_k^{\left(e\right)}}{\partial y^2}+\\ +2N_{12}^{N_{2p}}\left(u_{ij},v_{ij},w_i^{\left(e\right)},w_j^{\left(e\right)}\right)\frac{{\partial }^2{w}_k^{\left(e\right)}}{\partial x\partial y})w_r^{\left(e\right)}d\Omega ,\end{array}$ (33)

where

$\begin{array}{l}\left\{N^L\left(u_i^{\left(e\right)},v_i^{\left(e\right)}{w}_i^{\left(e\right)}\right)\right\}\left\{N_{11}^L,N_{22}^L,N_{12}^L\right\}=[A]\left\{{\epsilon }^{L_{0i}}\right\}+[B]\left\{{\chi }^i\right\},\\ \left\{{\epsilon }^{L_{0i}}\right\}={\left\{\frac{\partial u_i^{\left(e\right)}}{\partial x}-k_1{w}_i^{\left(e\right)},\frac{\partial v_i^{\left(e\right)}}{\partial y}-k_2{w}_i^{\left(e\right)},\frac{\partial u_i^{\left(e\right)}}{\partial y}+\frac{\partial v_i^{\left(e\right)}}{\partial x}\right\}}^T,\left\{{\chi }^i\right\}={\left\{-\frac{{\partial }^2{w}_i^{\left(e\right)}}{\partial x^2},-\frac{{\partial }^2{w}_i^{\left(e\right)}}{\partial y^2},-\frac{{\partial }^2{w}_i^{\left(e\right)}}{\partial x\partial y}\right\}}^T\\ \left\{N^{N_{2p}}\right\}=\left\{N_{11}^{N_{2p}},N_{22}^{N_{2p}},N_{12}^{N_{2p}}\right\}=\left[A\right]{\left\{{\epsilon }^{N_{2p}}\right\}}^T,\left\{M^{N_{2p}}\right\}=\left\{M_{11}^{N_{2p}},M_{12}^{N_{2p}},M_{12}^{N_{2p}}\right\}=\left[B\right]{\left\{{\epsilon }^{N_{2p}}\right\}}^T.\\ \left\{{\epsilon }^{N_{2p}}\left(u_{ij},v_{ij},w_i^{\left(e\right)},w_j^{\left(e\right)}\right)\right\}={\left\{{\epsilon }_{11}^{N_{2p}},{\epsilon }_{22}^{N_{2p}},{\epsilon }_{12}^{N_{2p}}\right\}}^T,{\epsilon }_{11}^{N_{2p}}=\frac{\partial u_{ij}}{\partial x}+\frac{1}{2}\frac{\partial w_i^{\left(e\right)}}{\partial x}\frac{\partial w_j^{\left(e\right)}}{\partial x},\\ {\epsilon }_{22}^{N_{2p}}=\frac{\partial v_{ij}}{\partial y}+\frac{1}{2}\frac{\partial w_i^{\left(e\right)}}{\partial y}\frac{\partial w_i^{\left(e\right)}}{\partial y},{\epsilon }_{12}^{N_{2p}}=\frac{\partial v_{ij}}{\partial x}+\frac{\partial u_{ij}}{\partial y}+\left(\frac{\partial w_i^{\left(e\right)}}{\partial x}\frac{\partial w_j^{\left(e\right)}}{\partial y}+\frac{\partial w_j^{\left(e\right)}}{\partial x}\frac{\partial w_i^{\left(e\right)}}{\partial y}\right).\end{array}$

Solution to the system of equations (31) can be found with the use of various approximation methods. Here, numerical implementation has been done with one mode. Thus, instead of the system of equations (31), we find the solution to one second-order ODE of the form

$y_1^{"}\left(t\right)+\alpha y_1\left(t\right)+\beta y_1^2\left(t\right)+\gamma y_1^3\left(t\right)=\mathrm{0,}$ (34)

where the coefficients α, β, γ are calculated using formulas (32)-(33), for i = j = k =r = 1. Thus, $\alpha ={\omega }_L^2,\beta ={\beta }_{11}^{\left(1\right)},\gamma ={\beta }_{111}^{\left(1\right)}$. The study of equations of the form (34) involved a lot of scientists [Mahmoud Bayat et al (2012, 2013); Iman Pakar et al. (2014, 2012); Mahdi Bayat et al. (2014), and others]. In this paper we have carry out numerical investigations by employing the classical Runge-Kutta method.

4 NUMERICAL RESULTS

In order to validate the results obtained by means of the proposed approach, a few of test problems are investigated first.

4.1 Results Validation

Task 1. The natural frequency of FG square shallow shells with movable simply supported edges and different values of the dimensionless parameter a/h = 10;5 is analyzed. Aluminum and Alumina FG mixture Al/Al ₂ O ₃ are considered as constituent materials. Material properties of the FG mixture used in the present study are taken as follows (see references [Reddy et al. (1999); Shen (2009)]):

$\begin{array}{l}Al:\begin{array}{c}\end{array}{E}_m=70Gpa,\begin{array}{c}\end{array}{v}_m=\mathrm{0.3,}\begin{array}{c}\end{array}{\rho }_m=2707kg/m^3;\\ A{l}_2{O}_3:\begin{array}{c}\end{array}{E}_c=380GPa,\begin{array}{c}\end{array}{v}_c=\mathrm{0.3,}\begin{array}{c}\end{array}{\rho }_c=3800kg/m^3.\end{array}$

The boundary conditions follow:

$\begin{array}{l}v=w=M_x={\psi }_y=N_x=0\begin{array}{c}\end{array}at\begin{array}{c}\end{array}x=\pm \frac{a}{2},\\ u=w=M_y={\psi }_x=N_y=0\begin{array}{c}\end{array}at\begin{array}{c}\end{array}y=\pm \frac{a}{2}.\end{array}$

To solve this problem, the following solution structure [Kurpa (2009)] is employed

$u=f_2{\Phi }_1,\begin{array}{cc}& \end{array}v=f_1{\Phi }_2,\begin{array}{cc}& \end{array}w=\omega {\Phi }_3,\begin{array}{cc}& \end{array}{\psi }_x=f_2{\Phi }_4,\begin{array}{cc}& \end{array}{\psi }_y=f_1{\Phi }_5,$

where

$\begin{array}{l}\omega (x,y)=\left(f_1{\wedge }_0{f}_2\right),\\ f_1=\frac{1}{2a}\left(a^2-x^2\right)\ge \mathrm{0,}\begin{array}{cc}& \end{array}{f}_2=\frac{1}{2b}\left(b^2-y^2\right)\ge 0.\end{array}$

Φ_i , i = 1, ..., 5 are undefined components represented as a truncated expansion of the complete system of functions:

${\Phi }_i={\displaystyle \sum _{k=1}^{k=N_i}{a}_k^{(i)}{\phi }_k^{(i)}}.$

In this paper, the power polynomials $\left\{{\phi }_k^{\left(i\right)}\right\}$ are chosen for the functions of such a system. In order to verify the accuracy of the presented results, convergence of the numerical solution is examined. As a result of the computational experiment, it has been found that convergence of the results (at least to the third decimal) occurs for the ninth-degree polynomials Φ₁, Φ₂, Φ₄, Φ₅, and for the tenth-degree polynomial Φ₃, and the results presented below were obtained for such a number of coordinate functions.

A comparison of the natural frequencies ${\Omega }_1={\lambda }_{1}h\sqrt{{\rho }_c/E_c}$ for various shallowness ratios and different thicknesses is carried out for $\frac{a}{R_y}=\frac{a}{R_x}=0$ (plate), $\frac{a}{R_y}=\frac{a}{R_x}$ (spherical shell), $\frac{a}{R_y}=0$ (cylindrical shell), and $\frac{a}{R_y}=-\frac{a}{R_x}$ (hyperbolic paraboloid).

The comparison of the obtained results for the side-to-thickness ratio a/h = 10 is shown in Table 1 and for a/h = 5 in Table 2.

Table 1
Comparison of the fundamental frequency parameters ${\Omega }_1={\lambda }_{1}h\sqrt{{\rho }_c/E_c}$ of square FG shallow shells with simply supported movable edges (Al/Al ₂ O ₃, a/h = 10).

$\frac{a}{R_y}$	$\frac{a}{R_x}$	k	RFM (CPT)	RFM (FSDT)	(CPT) [Alijani et al. (2011b)]	(FSDT) [Chorfi and Houmat (2010)]	(HSDT) [Matsunaga (2008)]
0	0	0	0.0597	0.0576	0.0597	0.0577	0.0578
		0.5	0.0505	0.0489	0.0506	0.0490	0.0492
		1	0.0455	0.0441	0.0456	0.0442	0.0443
		4	0.0395	0.0382	0.0396	0.0383	0.0381
		10	0.0380	0.0365	0.0380	0.0366	0.0364
0.5	0.5	0	0.0770	0.0753	0.0779	0.0762	0.0751
		0.5	0.0665	0.0652	0.0676	0.0664	0.0657
		1	0.0605	0.0593	0.0617	0.0607	0.0601
		4	0.0508	0.0496	0.0519	0.0509	0.0503
		10	0.0472	0.0462	0.0482	0.0471	0.0464
0	0.5	0	0.0642	0.0622	0.0648	0.0629	0.0622
		0.5	0.0546	0.0531	0.0553	0.0540	0.0535
		1	0.0494	0.0481	0.0501	0.0490	0.0485
		4	0.0423	0.0411	0.0430	0.0419	0.0413
		10	0.0403	0.0389	0.0408	0.0395	0.0390
0.5	-0.5	0	0.0582	0.0562	0.0597	0.0580	0.0563
		0.5	0.0493	0.0477	0.0506	0.0493	0.0479
		1	0.0444	0.0430	0.0456	0.0445	0.0432
		4	0.0385	0.0372	0.0396	0.0385	0.0372
		10	0.0370	0.0356	0.0380	0.0368	0.0355

Table 2
Comparison of the fundamental frequency parameters ${\Omega }_1={\lambda }_{1}h\sqrt{{\rho }_c/E_c}$ of the square FG shallow shells with simply supported movable edges (Al/Al ₂ O ₃, a/h = 5).

a/h = 5
b/Ry	a/Rx	Method	k = 0	k = 0,5	k = 1	k = 4	k = 10	k = ∞
0	0	RFM	0.211	0.180	0.162	0.139	0.132	0.108
		[Matsunaga (2008)]	0.212	0.182	0.164	0.138	0.131	0.108
0.5	0.5	RFM	0.2297	0.196	0.177	0.150	0.141	0.117
		[Matsunaga (2008)]	0.2301	0.200	0.182	0.151	0.142	0.117
1	1	RFM	0.275	0.237	0.215	0.177	0.164	0.140
		[Matsunaga (2008)]	0.274	0.243	0.223	0.186	0.169	0.139
0	0.5	RFM	0.214	0.183	0.165	0.141	0.133	0.109
		[Matsunaga (2008)]	0.215	0.186	0.168	0.141	0.133	0.110
0	1	RFM	0.223	0.191	0.173	0.146	0.137	0.114
		[Matsunaga (2008)]	0.224	0.194	0.177	0.148	0.138	0.114
-0.5	0.5	RFM	0.205	0.175	0.158	0.135	0.128	0.04
		[Matsunaga (2008)]	0.206	0.177	0.160	0.135	0.127	0.105
-1	1	RFM	0.191	0.163	0.148	0.126	0.119	0.097
		[Matsunaga (2008)]	0.192	0.165	0.149	0.125	0.118	0.098

The comparison shows that the results obtained using the refined first-order theory (RFM, FSDT) are almost the same as those reported in reference [Chorfi and Houmat (2010)]. A deviation from the results obtained by means of the theory of the higher order (HSDT) [Matsunaga (2008)] does not exceed 4%. Deviation results obtained by using the classical theory (RFM, CST) with the results of [Alijani et al. (2011a)] do not exceed 2%. In general, it should be noted that the classical theory, in most of cases, overestimates the fundamental frequencies compared with the refined theory.

4.2 Free Vibrations of Functionally Graded Shells of Complex Planforms

In order to present novel results and illustrate the versatility and efficiency of the proposed method, two free vibration problems are considered. Let us investigate a shallow shell with complex planform (see Fig. 1). The fixed geometrical parameters are as follows: h/2a = 0.1; b/2a ; a ₁/2a = 0.2; k ₁/k ₂ = (0;1, -1); b ₁/2a = (0.3; 0.35; 0.51); k ₁ = 2a/^_Rx ; k ₂ = 2a/^_Ry .

The properties of the FG mixture are the same as those presented in paper [Chorfi and Houmat (2010)], i.e.,

$\begin{array}{l}\left(\text{FG1}\right)Al/A{l}_2{O}_3:E_m/E_c=70/380GPa,v_m=v_c=\mathrm{0.3,}{\rho }_m/{\rho }_c=2707/3800kg/m^3;\\ \left(\text{FG2}\right)Al/Zr{O}_2:E_m/E_c=70/151GPa,v_m=v_c=\mathrm{0.3,}{\rho }_m/{\rho }_c=2707/3000kg/m^3.\end{array}$ (35)

Figure 1
Form and planform of the shells under consideration.

Suppose that the shell is clamped. Then, the solution structure may be taken in the following form

$u=\omega {\Phi }_1,\begin{array}{cc}& \end{array}v=\omega {\Phi }_2,\begin{array}{cc}& \end{array}w=\omega {\Phi }_3,\begin{array}{cc}& \end{array}{\psi }_x=\omega {\Phi }_4,\begin{array}{cc}& \end{array}{\psi }_y=\omega {\Phi }_5$ (36)

where ω = 0 is the equation of the border of the shell planform.

In order to realize the solution structure (36), one should construct the equation of the border ω = 0. Using the R-operations ∧₀, ∨₀ [Rvachev (1982)], the equation is built in the form

$\omega =\left(f_1{\wedge }_0{f}_2\right){\vee }_0{f}_3,$

where

$f_1=\left(\left(a_1^2-x^2\right)/2a_1\right)\ge 0$

is a vertical band bounded by straight lines x ± a ₁, and

$f_2=\left(\left(b_1^2-y^2\right)/2b_1\right)\ge 0$

is a horizontal band bounded by straight lines y ± b ₁. Finally,

$f_3=\left(1-\frac{x^2}{a^2}-\frac{y^2}{b^2}\right)\ge 0$

is a part of the plane within the ellipse.

Due to the doubly symmetric nature of the shell, numerical implementation is performed only for one-quarter of the investigated domain. Thus, the sequences of polynomials are chosen as follows

$\begin{array}{l}{\Phi }_1,{\Phi }_4:x,x^3,x{y}^2,x^5,x^3{y}^2,x{y}^4,x^7,x^5{y}^2,x^3{y}^4,x{y}^6,\cdots ;\\ {\Phi }_2,{\Phi }_5:y,x^{2}y,y^3,x^{4}y,x^2{y}^3,y^5,x^{6}y,x^4{y}^3,x^2{y}^5,y^7,\cdots ;\\ {\Phi }_3:\mathrm{1,}{x}^2,y^2,x^4,x^2{y}^2,y^4,x^6,x^4{y}^2,x^2{y}^4,y^6,\cdots .\end{array}$

Now, in order to investigate convergence of the natural frequencies, the computer experiments are carried out. It has been found that the third decimal is stabilized while maintaining the degree of approximating polynomial (11, 11, 14, 11, 11) which corresponds to the following number of coordinate functions for u, v, w, _^ψx , _^ψy: 21, 21, 36, 21, 21, respectively.

Fig. 2 shows the dependence of the natural frequencies ${\Omega }_L=4{\lambda }_1{a}^2\sqrt{{\rho }_c/E_c}$ versus the volume fraction exponent k. To check the reliability of the results, the calculation has been performed for the value of the parameter b ₁/2a = 0.51 In this case, the form shown in Fig. 1 is very close to elliptical. Results for this type of the shell have been compared to the similar form of the shell studied in reference [Chorfi and Houmat (2010)]. It is obvious that the results obtained for simple (ellipse, black and blue colors of curves) and complicated (dashed lines) forms coincide reasonably well. This fact allows us to validate our (numerically obtained) results. The same figure shows the effect of the volume fraction exponent k regarding the values of the natural frequency for spherical shells with radii of curvature equal to 2a/^_Rx = k ₁ = 0.5 and 2a/^_Ry = k ₂ = 0.5 The results have been obtained for two types of materials. As in the case of elliptical shells, the frequency values are substantially greater for the FG2 mixture than for the mixture of FG1.

Figure 2
Effect of volume fraction index k on the natural frequencies of the spherical and cylindrical panels.

The values of the natural frequency parameter ${\Omega }_L=4{\lambda }_1{a}^2\sqrt{{\rho }_c/E_c}$ for spherical, cylindrical, hyperbolic paraboloid shells and plates are presented in Table 3.

Table 3
Effect of volume fraction index k on the natural frequencies ${\Omega }_L=4{\lambda }_1{a}^2\sqrt{{\rho }_c/E_c}$ of the shallow shells with planform shown in Fig. 1 (a ₁/2a 0.2; b ₁/2a = 0.35).

$\frac{2a}{R_x}$	$\frac{2a}{R_y}$	FGM	k = 0	k = 1	k = 4	k = 10
0	0	FG1	0.8248	0.6457	0.5472	0.5109
		FG2	0.8248	0.7121	0.6633	0.6363
0.5	0.5	FG1	0.8707	0.6839	0.5760	0.5355
		FG2	0.8707	0.7518	0.6958	0.6692
0	0.5	FG1	0.8429	0.6604	0.5581	0.5203
		FG2	0.8429	0.7275	0.6747	0.6491
0.5	-0.5	FG1	0.8510	0.6681	0.5643	0.5254
		FG2	0.8510	0.7351	0.6813	0.6552

It follows from Table 3 those values of the natural frequencies decrease for all types of the shell curvature and material properties of the mixtures, provided that the volume fraction exponent k increases. Frequencies ‘asymptotically’ approach the frequencies of the metal shell or of the plate. It should be noted that for all values of k ∈ [0, 10] the spherical shells have the maximum value of the fundamental frequencies, while plates have the smallest value of the fundamental frequencies.

In order to carry out nonlinear analysis let us analyze the dependence between the amplitude and the ratio of nonlinear frequency to linear frequency. The backbone curves for clamped spherical shells made of FG1 material are shown in Fig. 3 while for those made of material FG2 are presented in Fig. 4. Geometrical parameters b ₁/2a and a ₁/2a are taken as follows: b ₁/2a = 0.35; a ₁/2a = 0.2.

Figure 3
Backbone curves for a clamped FG1 spherical shell (k ₁ = k ₂ 0.5; b ₁/2a = 0.35; a ₁/2a = 0.2).

Figure 4
Backbone curves for a clamped FG2 spherical shell (k ₁ = k ₂ 0.5; b ₁/2a = 0.35; a ₁/2a = 0.2).

Let us consider more shallow shells with curvatures k ₁ = 0.2, k ₂ = 0.2 and k ₁ = 0.2, k ₂ = 0, made of materials FG1 and FG2 provided that the remaining geometric parameters are the same as in the previous case. Effects of various values of volume fractions k on the ratio of nonlinear frequency to linear frequency and the amplitude for the spherical shell are shown in Table 4.

Table 4
Effect of volume fraction index k on the ratio ω_N /ω_L for clamped spherical shell (k ₁ = k ₂ 0.2; b ₁/2a = 0.35; a ₁/2a = 0.2).

$\frac{W_{\max }}{h}$	FGM	k = 0	k = 1	k = 4	k = 6	k = 10
0.25	FG1	1.0130	1.0121	1.0103	1.0083	1.0099
	FG2	1.0130	1.0125	1.0102	1.0119	1.0128
0.5	FG1	1.0504	1.0473	1.0419	1.0386	1.0398
	FG2	1.0504	1.0449	1.0430	1.0453	1.0455
0.75	FG1	1.1091	1.1057	1.0929	1.0864	1.0897
	FG2	1.1091	1.1059	1.0973	1.0991	1.1003
1	FG1	1.1868	1.1864	1.1639	1.1522	1.1558
	FG2	1.1868	1.1858	1.1679	1.1714	1.1727
1.25	FG1	1.2806	1.2831	1.2523	1.2380	1.2373
	FG2	1.2806	1.2771	1.2556	1.2587	1.2607
1.5	FG1	1.3873	1.3964	1.3562	1.3319	1.3323
	FG2	1.3873	1.3848	1.3559	1.3596	1.3609
1.75	FG1	1.5046	1.5216	1.4711	1.4417	1.4381
	FG2	1.5046	1.5037	1.4676	1.4708	1.4726
2	FG1	1.6333	1.6571	1.5959	1.5600	1.5519
	FG2	1.6333	1.6311	1.5887	1.5913	1.5924

Effects of various values of volume fractions k on ratio of nonlinear frequency to linear frequency and amplitude for the cylindrical shell with radii of curvature equal to 2a/_^Rx = k ₁ = 0.2 and 2a/_^Ry = k ₂ = 0, made of FG1 and FG2 materials, are shown in Table 5.

Table 5
Effect of volume fraction index k on the ratio ω_N /ω_L for clamped cylinder shell (k ₁ = 0.2; k ₂ = 0; b ₁/2a = 0.35; a ₁/2a = 0.2).

$\frac{W_{\max }}{h}$	FGM	k = 0	k = 1	k = 4	k = 6	k = 10
0.25	FG1	1.0143	1.0145	1.0113	1.0122	1.0119
	FG2	1.0143	1.0141	1.0132	1.0131	1.0131
0.5	FG1	1.0553	1.0563	1.0492	1.0472	1.0461
	FG2	1.0553	1.0546	1.0510	1.0507	1.0509
0.75	FG1	1.1197	1.1224	1.1075	1.1030	1.1004
	FG2	1.1197	1.1183	1.1106	1.1098	1.1103
1	FG1	1.2033	1.2096	1.1854	1.1773	1.1724
	FG2	1.2033	1.2017	1.1888	1.1873	1.1881
1.25	FG1	1.3029	1.3145	1.2800	1.2674	1.2594
	FG2	1.3029	1.3013	1.2825	1.2802	1.2812
1.5	FG1	1.4151	1.4336	1.3884	1.3707	1.3591
	FG2	1.4151	1.4141	1.3889	1.3856	1.3867
1.75	FG1	1.5372	1.5640	1.5079	1.4847	1.4689
	FG2	1.5372	1.5374	1.5055	1.5012	1.5022
2	FG1	1.6672	1.7032	1.6363	1.6073	1.5871
	FG2	1.6672	1.6689	1.6303	1.6247	1.6258

The related backbone curves for spherical and cylindrical shell, made of FG1 material are presented in Fig. 5. and Fig. 6, respectively.

Figure 5
Effect of volume fraction index k on the ratio ω_N /ω_L of spherical shell panel for FG1 material (k ₁ = 0.2; k ₂ = 0.2; b ₁/2a = 0.35; a ₁/2a = 0.2).

Figure 6
Effect of volume fraction index k on the ratio ω_N /ω_L of cylindrical shell panel for FG1 material (k ₁ = 0.2; k ₂ = 0; b ₁/2a = 0.35; a ₁/2a = 0.2).

A comparison of the backbone curves for cylindrical and spherical clamped panels for two kinds of materials is presented in Fig. 7 (volume fraction index k = 4).

Figure 7
Comparison of the backbone curves for the clamped cylindrical and spherical shells for materials FG1 and FG2 (b ₁/2a = 0.35; a ₁/2a = 0.2; kr1 = k ₁; kr2 = k ₂; k = 4)

Figure 8 shows the effect of the index k on the ratio of nonlinear to linear frequency. It can be seen that for materials FG1 this ratio increases to k = 1, and then it decreases to k = 6. However, if k > 6 then the ratio again increases. The same behavior is characteristic for the cylindrical shell made of material FG1 (see Table 5). But if the panels are made from materials FG2, then the ratio _^ωN /_^ωL increases to k = 6, and then it slowly begins to increase (see Tables 4, 5).

Figure 8
Effect of volume fraction index k on the ratio ω_N /ω_L of the spherical shell panel for material FG1 (k ₁ = 0.2; k ₂ = 0.2; b ₁/2a = 0.35; a ₁/2a = 0.2).

The carried out analysis of the backbone curves for the shells under consideration shows that these curves are sufficiently close to each other. However, the values of the nonlinear frequencies depend substantially on the volume fraction index k. Effects of volume fraction index k on the nonlinear frequencies of the clamped spherical and cylindrical shells are shown in Figure 9 and 10, respectively.

Figure 9
Effect of the volume fraction index k on nonlinear frequencies ${\Omega }_N={\omega }_N{a}^2\sqrt{{\rho }_c/E_c}$ of the spherical panel (FG1, k ₁ = k ₂ = 0.2; b ₁/2a =035; a ₁/2a = 0.2).

Figure 10
Effect of the volume fraction index k on nonlinear frequencies ${\Omega }_N={\omega }_N{a}^2\sqrt{{\rho }_c/E_c}$ of the cylindrical panel (FG1, k ₁ = 0.2; k ₂ = 0; b ₁/2a = 035; a ₁/2a = 0.2).

It should be emphasized that the panels with the highest values of nonlinear frequencies are made of pure ceramics.

The phase planes of the spherical and cylindrical considered shells for different initial conditions are shown in Figures 11 and 12, respectively.

Figure 11
Phase plane of the spherical shell for the volume fraction index k = 4 (k ₁ = k ₂ = 0.2; b ₁/2a = 035; a ₁/2a = 0.2) a) shell material-FG1; b) shell material FG2.

Figure 12
Phase plane of the cylindrical shell for the volume fraction index k = 4, (k ₁ = 0.2; k ₂ = 0; b ₁/2a = 035; a ₁/2a = 0.2) a) shell material-FG1; b) shell material FG2.

The proposed approach has allowed for a detailed analysis of the clamped cylindrical and spherical panels with the planforms shown in Fig. 1 being made of two kinds of FGM. It should be mentioned that owing to employment of the R-functions theory, we can easily pass from one geometric form to another and hence to investigate different boundary conditions of the shell using the same developed software. Our research shows that it is important that the desired solution to linear auxiliary tasks is presented in an analytical form. This is a significant factor that should be taken into account while using the presented approach to solve nonlinear problems.

5 CONCLUSIONS

This paper proposes a method of investigation of geometrically nonlinear free vibrations of functionally graded shallow shells of a complex planform. The method is based on the theory of the R-functions, Ritz variational method, Bubnov-Galerkin procedure, and Runge-Kutta method. The tests conducted for shells of square and elliptical planforms have proved the reliability and effectiveness of the presented method. Graphical and numerical results are obtained for shells having complicated shapes. In the future, the developed method is planned to be implemented into a multi-mode approximation.

References

• Alijani, F., Amabili, M., Karagiozis, K., Bakhtiari-Nejad, F. (2011b) Nonlinear vibrations of functionally graded doubly curved shallow shells. Journal of Sound and Vibration 330: 1432-1454.
• Alijani, F., Bakhtiari-Nejad, F., Amabili, M. (2011a). Nonlinear vibrations of FGM rectangular plates in thermal environments. Nonlinear Dynamic 66: 251-270.
• Amabili, M. (2008). Nonlinear Vibrations and Stability of Shells and Plates. Cambridge University Press, Cambridge.
• Awrejcewicz, J., Kurpa, L., Mazur, O. (2010). On the parametric vibrations and meshless discretization of orthotropic plates with complex shape. International Journal of Nonlinear Sciences and Numerical Simulation 11(5): 371-386.
• Awrejcewicz, J., Kurpa, L., Shmatko, T. (2013). Large amplitude free vibration of orthotropic shallow shells of complex shapes with variable thickness. Latin American Journal of Solids and Structure 10: 147-160.
• Awrejcewicz, J., Kurpa, L., Shmatko, T. (2015a). Investigating geometrically nonlinear vibrations of laminated shallow shells with layers of variable thickness via the R-functions theory. Composite Structures 125: 575-585.
• Awrejcewicz, J., Kurpa, L., Shmatko, T. (2015b). Vibration of functionally graded shallow shells with complex shape. In: Awrejcewicz, J., Kaźmierczak, M., Mrozowski, J., Olejnik, P. (Eds.) Dynamic Systems. Mathematical and Numerical Approaches 57-68.
• Bayat, M., Bayat, M., Pakar, I. (2014). Nonlinear vibration of an electrostatically actuated microbeam. Latin American Journal of Solids and Structures 11(3): 534-544.
• Bayat, M., Pakar, I. (2013). On the approximate analytical solution to non-linear oscillation systems. Shock and vibration 20(1): 43-52.
• Bayat, M., Pakar, I., Domairry, G. (2012). Recent developments of some asymptotic methods and their applications for nonlinear vibration equations in engineering problems: A review. Latin American Journal of Solids and Structures 9(2): 1-93.
• Bich, D.H., Dung, D.V., Nam, V.N. (2012). Nonlinear Dynamic analysis of eccentrically stiffened functionally graded cylindrical panels. Composite Structures 94: 2465-2473.
• Chorfi, S.M., Houmat, A. (2010). Non-linear free vibration of a functionally graded doubly-curved shallow shell of elliptical plan-form. Composite Structures 92: 2573-2581.
• Kurpa, L., Pilgun, G., Amabili, M. (2007). Nonlinear vibrations of shallow shells with complex boundary: R-functions method and experiments. Journal of Sound and Vibration 306: 580-600.
• Kurpa, L., Shmatko, T., Timchenko, G. (2010). Free vibration analysis of laminated shallow shells with complex shape using the R-functions method. Composite Structures 93: 225-233.
• Kurpa, L.V. (2009a). Nonlinear free vibrations of multilayer shallow shells with a symmetric structure and with a complicated form of the plan. J. Math. Sci. 162(1): 85-98.
• Kurpa, L.V. (2009b). R-functions Method for Solving Linear Problems of Bending and Vibrations of Shallow Shells. Kharkiv NTU Press, Kharkiv (in Russian).
• Kurpa, L.V., Mazur, O.S. (2010). Method of R-functions for investigating parametric vibrations of orthotropic plates with complex shape. Journal of Mathematical Sciences 171(4): 1-14.
• Kurpa, L.V., Mazur, O.S., Tkachenko, V.V. (2013). Dynamical stability and parametrical vibrations of the laminated plates with complex shape. Latin American Journal of Solids and Structures 10: 175-188.
• Kurpa, L.V., Mazur, O.S., Tkachenko, V.V. (2014). Parametric Vibrations of Multilayer Plates of Complex Shape. Journal of Mathematical Sciences 174(2): 101-114.
• Kurpa, L.V., Shmatko, T.V. (2014). Nonlinear vibrations of laminated shells with layers of variable thickness. Shell Structures: Theory and Applications 3: 305-308.
• Lam, K.Y., Li, H., Ng, T.Y., Chua, C.F. (2002). Generalized Differential quadrature method for the free vibration of truncated conical panels. Journal of Sound and Vibration 251: 329-348.
• Loy, C.T., Lam, K.Y., Reddy, J.N. (2008). Vibration of functionally graded cylindrical shells. Int. J. Mech. Sci. 41: 309-324.
• Matsunaga, H. (2008). Free vibration and stability of functionally graded shallow shells according to a 2D higher-order deformation theory. Composite Structures 84: 132-146.
• Neves, A.M.A., Ferreira, A.J.M., Carrera, E., Cinefra, M., Roque, C.M.C., Jorge, R.M.N., et al., (2013). Free vibration analysis of functionally graded shells by a higher-order shear deformation theory and radial basis functions collocation accounting for through the thickness deformations. Eur. J. Mech. A/Solids 37: 24-34.
• Pakar, I., Bayat, M. (2012). Analytical study on the non-linear vibration of Euler-Bernoulli beams. Journal of Vibroengineering 14(1): 216-224.
• Pakar, I., Cveticanin, L. (2014). Nonlinear free vibration of systems with inertia and static type cubic nonlinearities: an analytical approach. Mechanism and Machine Theory 77: 50-58.
• Pradyumna, S., Bandyopadhyay, J.N. (2008) Free vibration analysis of functionally graded curved panels using a higher-order finite element formulation. Journal of Sound and Vibration 318: 176-192.
• Reddy, J.N. (2009). Analysis of functionally graded plates. International Journal for numerical methods in engineering 47: 663-684.
• Reddy, J.N., Loy, C.T., Lam, K.Y. (1999). Vibration of functionally graded cylindrical shells. Int. J. Mech. Sci. 41: 309-324.
• Rvachev, V.L. (1982). The R-Functions Theory and Its Some Application. Naukova Dumka, Kiev (in Russian).
• Shen, H.S. (2009). Functionally Graded Materials of Plates and Shells. CRC Press, Florida.
• Sundararajan, N., Prakash, T., Ganapathi, M. (2005). Nonlinear free flexural vibrations of functionally graded rectangular and skew plates under thermal environments. Finite Elem. Anal. Des. 42: 152-168.
• Tornabene, F., Liverani, A., Caligiana, G. (2011). FGM and laminated doubly curved shells and panels of revolution with a free form meridian: a 2-D GDQ solution for free vibrations. Int. J. Mech. Sci. 53: 446-470.
• Tornabene, F., Viola, E. (2007). Vibration analysis of spherical structures elements using the GDQ method. Comput. Meth. Appl. 53: 1532-1560.
• Tornabene, F., Viola, E. (2009). Free vibration analysis of functionally graded panels and shells of revolution. Meccanica 44: 255-281.
• Woo, J., Meguid, S.A., Ong, L.S. (2006). Nonlinear free vibration of functionally graded plates. Journal of Sound and Vibration 289: 595-611.
• Xiang, X., Zheng, H., Guoyong, J. (2015). Free vibration of four-parameter functionally graded spherical and parabolic shells of revolution with arbitrary boundary conditions. Composites Part, B 77: 59-73.
• Zhao, X., Liew, K.M. (2009). Geometrically nonlinear analysis of functionally graded shells. Int. J. Mech. Sci. 51: 131-144.
• Zhu, S., Guoyong, J., Tiangui, Ye. (2014). Free vibration analysis of moderately thick functionally graded open shells with general boundary conditions. Composite Structures 117: 169-186.
• Zhu, S., Guoyong, J., Tiangui, Ye., Xueren, W. (2015) Free vibration analysis of laminated composite and functionally graded sector plates with general boundary conditions. Composite Structures 132: 720-736.

Publication Dates

• Publication in this collection
  Sept 2017

History

• Received
  09 Mar 2017
• Reviewed
  20 June 2017
• Accepted
  21 June 2017