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 strainstress 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 FourierRitz 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 (2009}a); ^{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 Rfunctions theory, variational Ritz method, BubnovGalerkin procedure, and RungeKutta 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 numericalanalytical method based on the Rfunctions 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
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 hightemperature 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
where P _{0}, P _{1} P _{1} P _{2} P _{3} 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
or equivalently,
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
According to the nonlinear firstorder shear deformation theory of shallow shells (FSDT), displacement components _{ ui } , u _{2}, u _{3} 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})]:
Strains ε = {ε_{11}; ε_{22}; ε_{12}}^{ T } at an arbitrary point of the shallow shell follow
where
and k _{1} = 1 / _{ Rx } , k _{2} = 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:
where
Strains χ = {χ _{11}; χ _{22}; χ _{12}}^{ T } are defined by the following formulas
where indicator δ is the tracing constant taking values 1 or 0 for the FSDT and CST, respectively.
The strain resultants N = (N _{11}, N _{22}, N _{12})^{ T } moment resultants M = (M _{11}, M _{22}, M _{12})^{ T } and shear stress resultants Q = (_{ Qx } , _{ Qy } )^{ T } are calculated by integration along Ozaxis, and they have the following forms
where
Transverse shear force resultants _{ Qx } , _{ Qy } are defined as
where
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 inplane force resultants N = (N _{11}, N _{22}, N _{12})^{ T } and moment resultants M = (M _{11}, M _{22}, M _{12})^{ T } in the framework of FSDT, taking into account the power law (1), are defined as follows
Mass density ρ is also estimated by integration along the shell thickness, to yield
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
where λ is a vibration frequency. Applying Hamilton’s principle, we get a variational equation in the following form
where U and T are strain and potential energy, respectively, and
whereas
Minimization of functional (19) is performed using the Ritz method, and the necessary sequence of coordinate functions is built with the help of the Rfunctions theory [^{Rvachev (1982}); ^{Kurpa (2009}b)].
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:
where linear differential operators _{ Lij } (i, j, = 1, 2, 3) are defined in the following way
The expressions of the nonlinear terms standing on the righthand side of the system (21)(23) are defined as follows
where:
The vector {_{ (N } } is defined by formula (8).
Note that the type of linear operators L
_{13}, L
_{23}, L
_{33} is simplified in the case of plates due to the condition k
_{1} = k
_{2} = 0, which implies
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
The functions _{ uij } , _{ vij } should satisfy the following system of differential equations
where
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:
where
The functions F _{1} and F _{2} are:
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
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 Rfunctions theory.
Substituting expressions (24)(25) for the functions u, v, w into equations of motion (21)(23) and applying the BubnovGalerkin procedure, we obtain the following system of nonlinear ordinary differential equations for the unknown functions _{ yr } (t), (r = 1, ..., n):
The formulas defining the coefficients
where
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 secondorder ODE of the form
where the coefficients α, β,
γ are calculated using formulas (32)(33), for i = j = k =r = 1. Thus,
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 _{2} O _{3} 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})]):
The boundary conditions follow:
To solve this problem, the following solution structure [^{Kurpa (2009})] is employed
where
Φ_{ i } , i = 1, ..., 5 are undefined components represented as a truncated expansion of the complete system of functions:
In this paper, the power polynomials
A comparison of the natural frequencies
The comparison of the obtained results for the sidetothickness ratio a/h = 10 is shown in Table 1 and for a/h = 5 in Table 2.


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 
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 firstorder 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. (2011}a)] 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 _{1}/2a = 0.2; k _{1}/k _{2} = (0;1, 1); b _{1}/2a = (0.3; 0.35; 0.51); k _{1} = 2a/_{ Rx } ; k _{2} = 2a/_{ Ry } .
The properties of the FG mixture are the same as those presented in paper [^{Chorfi and Houmat (2010})], i.e.,
Suppose that the shell is clamped. Then, the solution structure may be taken in the following form
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 Roperations ∧_{0}, ∨_{0} [^{Rvachev (1982})], the equation is built in the form
where
is a vertical band bounded by straight lines x ± a _{1}, and
is a horizontal band bounded by straight lines y ± b _{1}. Finally,
is a part of the plane within the ellipse.
Due to the doubly symmetric nature of the shell, numerical implementation is performed only for onequarter of the investigated domain. Thus, the sequences of polynomials are chosen as follows
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
The values of the natural frequency parameter


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 _{1}/2a and a _{1}/2a are taken as follows: b _{1}/2a = 0.35; a _{1}/2a = 0.2.
Let us consider more shallow shells with curvatures k _{1} = 0.2, k _{2} = 0.2 and k _{1} = 0.2, k _{2} = 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.

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 _{1} = 0.2 and 2a/_{ Ry } = k _{2} = 0, made of FG1 and FG2 materials, are shown in Table 5.

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.
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 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).
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.
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.
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 Rfunctions 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 Rfunctions, Ritz variational method, BubnovGalerkin procedure, and RungeKutta 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 multimode approximation.