1 INTRODUCTION
Functionally graded materials (FGM’s) are the new types of composite materials that have the unique feature of both high thermal and high mechanical properties. The constituent materials of FGMs (usually ceramic and metal) are combined such that the variation of their volume fractions is continuous, giving rise to smooth and gradual change of the mechanical and thermal properties. Hence in spite of the traditional composites that may suffer high amounts of thermal stresses at the interfaces, FGMs could withstand high temperatures while maintaining their structural integrities and providing desirable mechanical properties. These properties makes FGMs a suitable choice for 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 multi-mode 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 simply-supported, 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 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 two-dimensional 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 frequency-amplitude 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.
2 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 V_{m} and V_{c} represent the volume fractions of metal and ceramic, N is the volume fraction index and h is the thickness of the shell. The effective material properties of FG shell (P_{eff} ) including Young’s modulus,E, density, ρ and Poisson’s ratio, v, are defined as,
where P_{m} and P_{c} denote the material properties of metal and ceramic. Using (Eqs. 1) and (2) the effective material properties can be expressed as a function of z,
It is evident from (Eq. 3) that for N = 0, the volume fraction of ceramic become zero and so P_{eff} = P_{m} i.e. the material composition is pure metal and with the increase in N, the material varies from pure metal to pure ceramic.
2.1 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}),
Where
and,
where A, B and D are the sub-matrices of the stiffness matrix and are defined as,
where Q_{ij} 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 N_{s} and N_{q} are the normal stress resultant in axial,s, and circumferential directions,q, respectively. Also, N_{sθ} is the shear stress resultant.
2.2 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, W_{ij} (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 non-homogeneous differential equation in F that its particular solution, F_{p} , can be obtained as,
where f_{i} (t)’s are functions of W_{ij} (t)’s. The homogenous solution, F_{h} , can also be written as,
where , and are stress resultants that arise from the constraints on the edges of the shell and can be determined in terms of W_{i} (t)’s by applying the average in-plane boundary conditions 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., v = 0 and N_{s} = 0) as follows (^{Amabili et al., 1998}),
The continuity of v is also fulfilled via the following relation,
where ∂v / ∂θ 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., F = F_{p}+F_{h} where F_{p} and F_{h} are determined by (Eqs. 12) and (13). Then using (Eq. 11) , and 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,
where ω_{ij} s are linear natural frequencies corresponding to the modes used in (Eq. 11) and C_{ijmqs} and Q_{ijmq} are the constants that depend on the material and geometric properties of the cylindrical shell.
2.3 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, W_{i} (t)’s and their time derivatives, , in terms of W11 and . This pair, which is referred to as the master coordinates, corresponds to the fundamental mode of the shell. Using the definition, , the following fifth order polynomial expansion is used to approximate the slave coordinates:
where , and aijk and bijk 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 ū and . The final fifth order expressions that will be obtained for Wij(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:
where Sij, P1ij and P2ij 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) into (Eq. 19), the uncoupled modal ODEs associated with the fundamental NNM of the cylindrical shell can be obtained as follows:
where,
with ki s being constant terms dependent on the material and geometric properties of the shell.
2.4 Solution by HAM
An efficient and accurate solution can be provided for (Eq. 22) using the HAM (^{Liao, 2003}). For this purpose, the transformation, τ = ωNLt and X(τ) = W1(t) - δ with ω and δ 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 and , it is reasonable to choose the following function as an initial guess of the solution,
The nonlinear operator, N is then defined as,
where q is the embedding parameter that varies from 0 to 1 and Ω(q), ∆(q) and Φ(τ, q) are the functions that represent the solution of the ODE when q = 1. The following linear operator is also defined in such a way that it’s homogeneous solution is cos(τ)
with ω0 being the first order approximation of ωNL. The zeroth-order deformation equation can now be defined as,
where K is the auxiliary undetermined constant parameter that if chosen properly may improve the convergence of the solution. Next, Ω(q), ∆(q) and Φ(t, q) are expanded in the Taylor series as,
which converges to the exact solutions of ωNL, δ and X(τ) when q = 1. That is to say,
where Xm(τ), δm and ωm are the unknowns that could be defined via constructing the higher-order deformation equations by differentiating the zeroth-order deformation equation m times with respect to q, then dividing it by m! and finally setting q = 0. The resulted m’th order deformation equations is,
where Xm = 0 for m = 1 and Xm = 1 for m > 1 and,
(Eq. 32) is a linear differential equation that along with the initial conditions, Xm(0) = 0 and 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(τ)), leading to algebraic equations in ω_{m} and δ_{m} . In numerical analysis of the present problem, it is found that the values of δ_{m} 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. Hence δ_{m} s are set to zero, which leads to the following relation for ω _{0}:
where . The first order approximation of the displacement is also obtained as:
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.
3 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, a1 defined in (^{Mahmoudkhani et al., 2011}) by , is calculated for a cylindrical shell made up of the stainless steel and the silicon nitride with n1 = 5,L/R = 2,h/R = 0.01,R = 0,01m, and different values of N. The results are compared with the result of (^{Mahmoudkhani et al., 2011}), as shown in Figure 2. A close agreement can be seen between the results which can ensure the validity of the present solution.
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 n1 = 10,L/R = 2, h/R = 0.004,R = 0.5m. The shell is also assumed to be composed of stainless steel and silicon nitride with the material properties given in Table 1. For the initial conditions (ICs), W11(0) = a11 and are used and the remaining ICs, required for solution of (Eqs. 19), are obtained using the IM approach to determine W_{ij} (0) ((i, j) ≠ (1, 1)) in terms of W _{11} and (i.e., the fifth order version of (Eq. 21)). The results, given in Figure 3 match well for a _{11}/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 a _{11} = h grows to 2.
Material | P _{0} | P _{-1} | P _{1} | P _{2} | P _{3} | |
---|---|---|---|---|---|---|
Stainless steel | E | 201.04e+9 | 0 | 3.079e-4 | -6.534e-7 | 0 |
α | 12.330e-6 | 0 | 8.086e-4 | 0 | 0 | |
ν | 0.3262 | 0 | -2.002e-4 | 3.797e-7 | 0 | |
ρ | 8166 | 0 | 0 | 0 | 0 | |
Silicon nitride | E | 348.43e9 | 0 | -3.07e-4 | 2.160e-7 | -8.946e-11 |
α | 5.8723e-6 | 0 | 9.095e-4 | 0 | 0 | |
ν | 0.2400 | 0 | 0 | 0 | 0 | |
ρ | 2370 | 0 | 0 | 0 | 0 |
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 k = -1 might be proper choice for the analysis. Time histories obtained by different order HAM for a _{11}/h = 2 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 a _{11}/h = 6 are also given in Figure 5-(b), which 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 index, N. The variation of the nonlinear frequencies (ω_{NL} ) with nondimensional amplitude (a _{11}/h) is given in Figure 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, ω_{NL} /ω _{11} v.s. a _{11}/h is depicted in Figure 6-(b). The slight increase in the bending of the curves is greater for intermediate values of N (N = 1,5). Similar study is also performed in Figure 7 for other system parameters taken as, n _{1}.= 8,L/R = 2,h/R = 0.0001, R =1m The nonlinearity is mostly hardening in this case, which is completely different from the previous one. The curves obtained in Figure 7-(b), also show the considerable effect of N on the nonlinearity. As is evident, for some intermediate values of N (N = 5,10), ω_{NL} / ω _{11} initially decrease with amplitude and then, after a certain amplitude, begins to rise. For other values of N, the most strong nonlinearity (i.e., bending of the backbone curve) occurs for N = 0 which decreases as N grows.
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. These surfaces depict the variation of the linear modal amplitudes, W _{01}, W _{03} and W _{21} with W _{11} and .
The results presented in Figure 8 for the first case study, show that W _{01} and W _{03} have a small dependence on although they considerably varies with W _{11}. It is in contrast with W _{21}, which has comparable dependency on both W _{11} and . However, the contribution of the mode associated with W21 to the NNM is not considerable as the magnitude of W21 is much smaller than the magnitudes of W01 and W03. 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, is seen to have more influence, in this case, on W _{01} and W _{03}, although the effect of W _{11} is still dominant. Moreover, W _{21} 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 on W _{21}, and also considering that W _{11} and W _{21} 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, θ is set to zero and for the circumferential direction, s is set to L/2. The linear mode shapes are also included in these figures to illustrate the effects of the nonlinearity on the mode shapes. Results in Figure 10 correspond to the first case study, whose corresponding invariant surfaces are presented in Figure 8. The mode shapes at t = 0 where is zero, are shown in Figure 10-(a). Slight differences between the linear and nonlinear mode shapes, both for longitudinal and circumferential directions, can be recognized, although the overall shapes are the same. Mode shapes at t = 0.25T_{NL} where T_{NL} = 2π / ω_{NL} is depicted in Figure 10-(b). The displacement amplitude (W _{11}) is nearly zero at this point, while is close to its maximum value. The nonlinear mode’s displacement distribution in both longitudinal and circumferential direction deviate considerably from their linear counterparts such that the number of circumferential half-wave numbers are doubled at this instant. This happens due to the fact that the amplitude of W _{21}, which corresponds to the mode with 2n _{1} circumferential half-wave number, becomes dominant at this time. However, since W _{21} have much lower amplitudes in most times than the modes with n _{1} circumferential half-wave numbers (see Figure 8), thus the nonlinear mode shape seen in Figure 10-(b) can be observed only in a small interval of time. The mode shapes at t = 0.5T_{NL} are also similar to those at t = 0, especially in the circumferential direction.
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 t = 0 and 0.5T_{NL} the mode shapes are similar to the trapezoidal form rather than the sinusoidal form. This is in fact due to larger amplitudes of W _{21} as is seen in Figure 9. Moreover, the circumferential wave numbers at t = 0.25T_{NL} (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 2n _{1} circumferential half-wave number is much more pronounced.
4 CONCLUSIONS
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 frequency-amplitude 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.