Nonlinear Vibration and Mode Shapes of FG Cylindrical Shells

The nonlinear vibration and normal mode shapes of FG cylindrical shells are investigated using an efficient analytical method. The equations of motion of the shell are based on the Donnell’s nonlinear shallow-shell, and the material is assumed to be gradually changed across the thickness according to the simple power law. The solution is provided by first discretizing the equations of motion using the multi-mode Galerkin’s method. The nonlinear normal mode of the system is then extracted using the invariant manifold approach and employed to decouple the discretized equations. The homotopy analysis method is finally used to determine the nonlinear frequency. Numerical results are presented for the backbone curves of FG cylindrical shells, nonlinear mode shapes and also the nonlinear invariant modal surfaces. The volume fraction index and the geometric properties of the shell are found to be effective on the type of nonlinear behavior and also the nonlinear mode shapes of the shell. The circumferential half-wave numbers of the nonlinear mode shapes are found to change with time especially in a thinner cylinder.

Latin American Journal of Solids and Structures 14 (2017) 422-440 use in aerospace structures such as outer skins of high speed vehicles, which are commonly in the form of cylindrical shells.The analysis of functionally graded (FG) cylindrical shells is therefore of critical importance and could lead to more efficient structures in the future.Vibration analyses are especially one of the main and critical considerations and thus should be given a high priority in designing cylindrical shell used in aerospace application.Numerous studies have been performed so far on the linear vibration of FG cylindrical shell considering the influence of different parameters such as the edge conditions, temperature rise and follower axial forces (Loy et al. 1999, Pradhan et al. 2000, Haddadpour et al. 2007, Torki et al. 2014).Vibration amplitudes of cylindrical shell may not, however, small enough to be accurately predicted by the linear analysis and thus the geometric nonlinearity should be accounted for.This has been the main concern of many studies in recent years.The work of Mahmoudkhani et al. (2011) is one of the first attempt in this area which is devoted to the analysis of nonlinear primary resonance response of FG cylindrical shells.The multimode Galerkin's procedure together with the method of multiple scales (MMS) is used to obtain the frequency-response curves and identify the ranges of excitation amplitudes and frequencies in which the multi-mode response occurs.The effects of temperature and axial loads on the free and forced nonlinear vibrations of FG cylindrical shells is investigated by Sheng and Wang (2013) and Bich and Xuan Nguyen (2012), using again the Galerkin's method.Rafiee et al. (2013) investigated the piezoelectric FG cylindrical shell under aerodynamic and thermal loads.Sofiyev (2016) studied the vibration of orthotropic FG cylindrical shells, giving an expression by the Jacobian elliptic function for the nonlinear frequencies.The homotopy perturbation method (HPM) is used to solve the single time equation obtained from the Galerkin's method.Cylindrical shells with different simplysupported, clamped and free boundary conditions are studied by Strozzi and Pellicano (2013).They used the Chebyshev orthogonal polynomial for discretization of the equations.Jafari et al. (2014) studied the FG cylindrical shell embedded with a piezoelectric layer and used the Multi-mode Galerkin and Runge-Kutta method to provide the solution.More detailed information and references on this subject may be found in the review paper by Alijani and Amabili (2014).
The method of solution adopted in most studies as those mentioned above, begins with discretizing the governing partial differential equations (PDEs) with the modal expansion method and then uses either the analytical or numerical methods to determine the periodic response.Among the analytical methods, the MMS is most frequently used due to its established routine.Lower order MMS may not, however, accurate enough as the oscillation amplitude increases.Using higher order MMS may also require heavy mathematical manipulation.Some other analytical methods based on the homogony method such as HPM or Homotopy analysis method (HAM) can give more accurate first order approximation (He 1999, Liao 2003) and their extension to higher orders can be readily performed.Hence they are proper choices for accurate and efficient solution of time equations especially when applied to a single equation.The single time equation obtained by using a single linear mode in the Galerkin's method is not, however, reliable and may yield erroneous result.Alternatively, the nonlinear normal modes (NNM) can be employed to decouple the set of nonlinear equations obtained by the multi-modal Galerkin's method.The single nonlinear equation that will be obtained in this way, could correctly describe the nonlinear behavior of the system.
The concept of NNM is first introduced by Rosenberg (1962) as a periodic synchronous oscillation of all degrees of freedom (DOF) of a nonlinear system.Since then, extensive studies have been Latin American Journal of Solids and Structures 14 (2017) 422-440 performed both on their applications and also their methods of computation (Vakakis 1997, Pierre et al. 2006, Kerschen et al. 2009, Peeters et al. 2009).The invariant manifold (IM) approach, proposed by Shaw and Pierre (1993) is one of these methods that defines the NNM as a twodimensional invariant manifold in phase space.The method extends the invariance property of the linear normal modes to nonlinear systems and may be considered as the generalization of the Rosenberg's definition (Peeters 2010).It is to be noted that in contrast to the linear normal modes, the NNMs are not orthogonal, but due to their invariance property they could still be suitable for accurate and efficient order reduction in nonlinear oscillatory systems (Kerschen et al. 2009).Moreover the mode shapes corresponding to NNMs would change with amplitude and thus with time (Touze et al. 2004).
IM method is applied in the present study to compute the NNMs of the FG cylindrical shell.For this purpose, the equations of motion, which are based on the Donnell's shell theory, are first discretized using the multi-mode Galerkin's method.The obtained ordinary differential equations are then decoupled using the IM method and then solved by HAM to determine the frequencyamplitude relation.Numerical results are obtained for backbone curves and nonlinear invariant surfaces.Nonlinear mode shapes are also presented by depicting the out-of-plane displacement variations in axial and circumferential directions.

FORMULATION
According to the simple power low, the variation of volume fractions of the ceramic and metal through the thickness direction ( z ) can be described as, where m V and c V represent the volume fractions of metal and ceramic, N is the volume fraction index and h is the thickness of the shell.
It is evident from Eq. (3) that for 0 N = , the volume fraction of ceramic become zero and so eff m P P = i.e. the material composition is pure metal and with the increase in N , the material varies from pure metal to pure ceramic.

Governing Equations
Considering the coordinate system shown in Figure 1 , the equation of motion and the compatibility equation of the FG cylindrical shell can be expressed in terms of redial displacement, w and the airy stress function, F as (Jansen 2008), and, where A , B and D are the sub-matrices of the stiffness matrix and are defined as, where ij Q are the components of the 3 3 ´ stiffness matrix based on the plane-stress Hook's law whose nonzero components are defined as, It must be noted that to drive the compatibility equation, Eq.( 5), the in-plane inertial terms should be neglected and the stress function, F , would be defined in the way to satisfy the in-plane equations of motion as, In which s N and N q are the normal stress resultant in axial, s , and circumferential directions, q , respectively.Also, s N q is the shear stress resultant.

Discretization of PDEs
The nonlinear PDEs obtained for the equations of motion (i.e., Eqs. ( 4) and ( 5)) are first discretized using the Galerkin's method with appropriate trial functions.Various approximations are used for this purpose in the literatures, among them the function used by (Amabili et al. 1998) contains the minimum number of the linear modal functions needed for accurately predicting the softening behavior of the shell.This function is defined as, where, ( ) ij W t 's are the unknown generalized coordinates to be determined.Substituting the displacement function, Eq. ( 11), into the compatibility equation, Eq. ( 5), yields a linear nonhomogeneous differential equation in F that its particular solution, p F , can be obtained as, where ( ) The homogenous solution, h F , can also be written as,

R z,w ,v z,w s,u h L
Latin American Journal of Solids and Structures 14 (2017) 422-440 where s N , N q and s N q are stress resultants that arise from the constraints on the edges of the shell and can be determined in terms of ( ) i W t 's by applying the average in-plane boundary condi- tions and the continuity requirement of circumferential displacement, v (Amabili et al. 1998).The above equation is not a general homogeneous solution and is chosen to satisfy the average in-plane boundary conditions which are defined for the classic simply supported cylindrical shell (i.e., and 0 s N = ) as follows (Amabili et al. 1998), The continuity of v is also fulfilled via the following relation, where v q ¶ ¶ can be related to w and F as, In the above equation F must be replaced by the complete solution of compatibility equation, i.e., where p F and h F are determined by Eqs. ( 12) and (13).Then using Eq. ( 11) s N , N q and s N q can be obtained in terms of the generalized coordinates as, Finally, introducing Eqs. ( 11), ( 12) and ( 13) into Eq.( 4), and applying the Galerkin's method, four nonlinear ODEs will be obtained as, Latin American Journal of Solids and Structures 14 (2017) 422-440 where ij w s are linear natural frequencies corresponding to the modes used in Eq. ( 11) and ijmqs C and ijmq Q are the constants that depend on the material and geometric properties of the cylindrical shell .

Extraction of Nonlinear Normal Modes
Eqs. ( 19) include coupled quadratic and cubic nonlinear terms that may be decupled via the corresponding NNMs.To extract the NNMs, the IM approach is used by expressing all the generalized coordinates, ( ) , the following fifth order polynomial expansion is used to approximate the slave coordinates: ( 15) 0 0 , ( , ) (0,1),(0, 3),(1,1),(2,1), where , and ijk a and ijk b are the unknown constants that could be determined by substitution of Eq. ( 20) into the state-space form of the Eqs.( 19) and equating the coefficients of like powers of u and v .The final fifth order expressions that will be obtained for ( ) ij W t are lengthy and may not be given here.But to give an insight to the characteristics of the NNMs, the third order approximation is presented as follows:  2) 2 (4 ) ( 4) ( , ) (0,1),(0, 3),(2,1).
where ij S , 1 ij P and 2 ij P are constant terms that depend on the system parameter values.It is evident from Eq. ( 21) that the fundamental NNM of the cylindrical shells may depend on both amplitude and velocity.This behavior, which is the common characteristic of nonlinear systems with quadratic nonlinearities, implies that every point of the structure in a certain NNM does not reach the zero-deformation state at the same instant (Nayfeh and Nayfeh 1995).Substituting Eq. ( 21) with i k s being constant terms dependent on the material and geometric properties of the shell.

Solution by HAM
An efficient and accurate solution can be provided for Eq. ( 22) using the HAM (Liao 2003).For this purpose, the transformation, -with w and d being the nonlinear frequency and the constant drift amplitude respectively, are used as follows, The unknown drift is included in the response, since Eq. ( 22) has quadratic and quartic terms.Considering the initial condition of (0) X W = and (0) 0 X =  , it is reasonable to choose the following function as an initial guess of the solution, 0 ( ) cos( ).
The nonlinear operator, Ν is then defined as, where q is the embedding parameter that varies from 0 to 1 and ( ) q W , ( ) q D and ( ; ) q f t are the functions that represent the solution of the ODE when 1 q = .The following linear operator is also defined in such a way that it's homogeneous solution is cos( ) t , with 0 w being the first order approximation of NL w .The zeroth-order deformation equation can now be defined as, Latin American Journal of Solids and Structures 14 (2017) 422-440 where k is the auxiliary undetermined constant parameter that if chosen properly may improve the convergence of the solution.Next, ( ) q W , ( ) q D and ( ; ) t q f are expanded in the Taylor series as, which converges to the exact solutions of NL w , d and ( ) X t when 1 q = .That is to say, (1) where ( ) m X t , m d and m w are the unknowns that could be defined via constructing the higher-order deformation equations by differentiating the zeroth-order deformation equation m times with re- spect to q , then dividing it by !m and finally setting 0 q = .The resulted m'th order deformation equations is, where Eq. ( 32) is a linear differential equation that along with the initial conditions, (0) 0 could be readily solved to obtain higher order approximations of the solution.This is performed along with eliminating the secular terms (i.e., coefficients of cos( ) t ), leading to algebraic equations in m w and m d .In numerical analysis of the present problem, it is found that the values of m d are negligible for various system parameters.This is because that the coefficients of quadratic and quartic terms in Eq. ( 22) are much smaller than those of cubic terms and can be neglected.
The first order solution presented by Eqs. ( 33) and ( 34), is adequately accurate for the range of initial displacement amplitudes considered in the present study.Arbitrarily higher order solutions can however be readily obtained using the mentioned procedure, albeit their corresponding expressions are so lengthy to be given here.

NUMERICAL RESULTS AND DISCUSSION
Numerical results are presented in this section to examine the accuracy of the IM method used for extraction of NNMs and also the solution of the obtained single nonlinear equation by the HAM.The nonlinear frequency-amplitude of FG cylindrical shells and also the nonlinear mode shapes of the cylindrical shell are also presented and compared with their linear counterparts.
For the verification process it should considered that the solution provided in the present study contains three types of approximation.The first approximation corresponds to the discretization of the partial differential equations (i.e., Eqs. ( 4) and ( 5)) which have led to the discretized equations defined by Eq. ( 19).The second approximation corresponds to the application of the IM approach to Eq. ( 19), and the third one is related to the analytical solution provided by using the HAM.Therefore, to assess the validity or accuracy of the results obtained by each of these approximations, three numerical studies are performed in Figures 2, 3 and 5. To ensure the validity and convergence of the discretization process used for obtaining Eq. ( 19), the nonlinearity index, 1 a defined in Mahmoudkhani et al. (2011) by , is calculated for a cylindrical shell made up of the stainless steel and the silicon nitride with 1 5 / 2, / 0.01, 0.1m , Next the accuracy of the NNM extracted by the IM approach from the discretized equations is examined in Figure 3, by comparing the fundamental nonlinear mode's response, obtained by the numerical solution of Eq.( 22), with the results of numerical simulation of the set of four equations given in Eq. ( 19).The system parameters used for simulations are The shell is also assumed to be composed of stainless steel and silicon nitride with the material properties given in Table 1.(i.e., the fifth order version of Eq. ( 21)).The results, given in Figure 3 match well for 11 / a h lower than 1.5 even when the third order IM is used.For higher values of initial displacement, it is seen that the fifth order IM remains accurate when 11 / a h grows to 2.  Reddy and Chin (1998).( 22)) by the numerical solution of set of four equations, Eq. ( 19), for (a) The accuracy of the HAM solution is also assessed in Figure 5 by comparing the numerical solution of Eq. ( 22) with the HAM solution.The system parameters are the same as those used in the previous study.To improve the convergence of the solution, the variations of the estimated nonlinear frequency and also the response amplitude at an arbitrary time instant with k are depicted first in Figure 4 Different-order HAM approximations are included in this figure.The best accuracy may be obtained by choosing k from the region where the variation of the solution with k is minimal.Hence shows that the fourth-order HAM have completely matched with the numerical result.The HAM, then is used to calculate the nonlinear backbone curves of the cylindrical shell with the same geometric properties of the previous study and different values of the volume fraction in-  6-(a).As N increases it is seen that the corresponding frequencies increase due to higher contribution of the ceramic part in the FGM.All the curves also show the softening behavior.
To better illustrate the effect of N on the strength of the softening nonlinearity, ).Similar study is also performed in Figure 7 for other system parameters taken as, 1 8 / 2, / 0.001, 1m , The nonlinearity is mostly hardening in this case, which is completely different from the previous one.The curves obtained in Figure 7  , The invariant nonlinear modal surfaces corresponding to the fundamental NNM of the shell are next obtained in Figures 8 and 9, for the two cases of the cylindrical shells studied above.W .This implies that the axisymmetric modes are mainly responsible for the softening behavior observed in Figure 6.The invariant surfaces of the second case study are shown in Figure 9. Comparing with the result of Figure 8, 11 W  is seen to have more influence, in this case, on 01 W and 03 W , although the effect of 11 W is still dominant.Moreover, 21 W has the greatest amplitude and thus the corresponding mode is the most dominant mode in the NNM compared to the axisymmetric modes.This mode is, therefore, the main cause of the hardening type nonlinearity seen in Figure 7. Also, due to the considerable effect of 11 W  on 21 W , and also considering that 11 W and 21 W correspond to the modes with one longitudinal half wave number and different circumferential wave numbers, it may be inferred that the nonlinear mode shape may significantly vary with time along the circumferential direction.This is examined in the next numerical studies in Figures 10 and 11.In these figures, the nonlinear mode shapes of the cylindrical shell is compared with their linear counterparts.
As is previously mentioned, the nonlinear mode shapes may not remain the same at different instants.Hence the mode shapes in Figures 10 and 11 are presented for three different time instants.Moreover at each time instants, the displacement distributions in longitudinal and circumferential directions are given in the same two dimensional plot.For the longitudinal distribution, q is set to zero and for the circumferential direction, s is set to / 2

L
. The linear mode shapes are also includ- , Mode shapes of the second case study, shown in Figure 11, exhibits a rather different nonlinear mode shapes in comparison to the first case.The displacement distribution in the longitudinal direction is nearly the same as the linear mode shapes, at different instants of time.However, the distributions in the circumferential direction are very different from the linear modes, such that at 0 t = and 0.5 NL T the mode shapes are similar to the trapezoidal form rather than the sinusoidal form.This is in fact due to larger amplitudes of 21 W as is seen in Figure 9.Moreover, the circumfer- ential wave numbers at 0.25 NL t T = (where the velocity has its maximum value), is twice the wave number of the linear mode shape.This is similar to what is observed for the first case in Figure 10-(b), except that the mode's maximum amplitude is much larger in this case and the nonlinear mode shape with 2 1 n circumferential half-wave number is much more pronounced.The nonlinear mode shapes of the FG cylindrical shell was calculated using the fifth-order IM approach and was employed along with the HAM to provide an expression for the frequencyamplitude relation of FG cylindrical shells.For this purpose, the governing PDEs were discretized using the most influential linear modes of the structure.These equations were then decoupled using the NNM, leading to a single nonlinear ordinary equation for the fundamental mode.The HAM was applied to solve this equation which was shown to give an accurate result even with the first approximation, in the reasonable range of displacement amplitudes considered in the present study.Numerical studies were performed to examine the accuracy of the IM and HAM.Moreover the backbone curves, the nonlinear invariant surfaces and also the nonlinear mode shapes of the FG cylindrical shell were obtained for two different thickness to radius ratios.The softening nonlinearity with minimal dependence on N was observed for the thicker shell.The thinner cylinder, however, exhibited the hardening behavior with the most nonlinearity occurred for lower values of N .It was found that the axisymmetric mode with twice circumferential wave number (i.e., mode (2,1)) is the dominant mode, compared to the axisymmetric modes, that contributes to the NNM in the thinner cylinder.Moreover, the dependence of (2,1) mode's amplitude on both the displacement and velocity of the mode (1,1) was found to be comparable, in contrast to the axisymmetric modes that were mostly dependent on the displacement.Hence in times that the displacement and velocity of the (1,1) mode reach their zero and maximum values, respectively, the amplitude of the mode (2,1) becomes dominant leading to the nonlinear mode shapes with circumferential wave number twice those observed for the linear mode shapes.This was in fact more pronounced for the thinner cylinder.

Figure 3 :
Figure 3: Comparison of the numerical solution of the single equation obtained by the third and fifth order IM (Eq.

1k
= -might be proper choice for the analysis.Time histories obtained by different order HAM for 11 / 2 a h = is compared in Figure 5-(a) with the numerical solution.Excellent agreement is seen in Figure 5-(a) with even first order HAM.To emphasize on the accuracy of the HAM for higher initial displacement, results obtained for 11 / 6 a h = are also given in Figure 5-(b), which

Figure 4 :
Figure 4: Variation of the frequency ratio and the response amplitude at time 0 0.004 t = sec with the auxiliary parameter k .
in Figure 6-(b).The slight increase in the bending of the curves is greater for intermediate values of N ( 1, 5 N =

Figure 8 :
Figure 8: The nonlinear invariant modal surface obtained for fundamental NNM by fifth order IM (a) 01 / W h (b) 03 / W h (c) 21 / W h v.s.11 / W h and 11 / NL W h w

Figure 9 :Figure 10 :
Figure 9: The nonlinear invariant modal surface obtained for fundamental NNM by fifth order IM 01 / W h (b) 03 / W h (c) 21 / W h v.s.11 / W h and 11 / NL W h w

Figure 11 :
Figure 11: Fundamental nonlinear mode shape of the cylindrical shell in different instants for 1 s and their time derivatives, ( )

Latin American Journal of Solids and Structures 14 (2017) 422-440 into
Eq. (19), the uncoupled modal ODEs associated with the fundamental NNM of the cylindrical shell can be obtained as follows:

Table 1 :
Material properties of FGM's from N .The variation of the nonlinear frequencies ( NL w ) with nondimensional amplitude ( 11 / Latin American Journal of Solids and Structures 14 (2017) 422-440 dex,