Nonlinear Parametric Dynamics of Bidirectional Functionally Graded Beams

In this paper, the parametric dynamics of bidirectional functionally graded (BDFG) beams subjected to a time-dependent axial force are studied. *e material properties of beam which vary along both thickness and axial directions follow the power law, and four different distribution patterns are considered. *e coupled nonlinear partial differential equations describing the longitudinal-transverse displacements and the shear deformation are derived using Hamilton’s principle based on Timoshenko beam theory. *e Galerkin scheme is employed to discrete the continuous model resulting in a multiple degree-of-freedom system, namely, the reduced order model. *e nonlinear parametric response of the beam is obtained by solving the discrete system numerically, and the frequencyand force-response curves are constructed by tracing the period motion using the pseudoarclength continuation technique. Numerical results are presented to examine the effects of system parameters, e.g., gradient parameters, magnitude and frequency of external excitation, and damping coefficients. Cyclic-fold bifurcation and branch points of the period motion are spotted in parametric resonance of the BDFG beam. Results show that the asymmetrical material distribution in thickness direction of beam leads to the asymmetry of dynamic responses. Moreover, the gradient of material in axial direction has more significant effect on the dynamic features of BDFG beam than that in the thickness direction.


Introduction
With the rapid development of engineering construction, various composite materials have been increasingly used such as the classic laminated composite and sandwich structures [1][2][3][4][5][6][7] which have advantages of high stiffness and strength-to-weight ratios and long fatigue life. However, classic laminated composites and sandwich structure may produce great interlaminar stresses between layers, when it is deformed. Functionally graded materials (FGMs) are a class of composites which possess continuously varying microstructure and mechanical properties without internal boundaries and hence can moderate the interfacial stress concentrations in classic laminated composites and sandwich structure. e gradual variation which results in a very efficient material can be designed for the specified purpose of achieving thermal resistance or corrosion resistance, etc. In the last few decades, the conventional unidirectional FGMs in which the gradual variation of material only exists in one direction [8][9][10][11] have been already widely reported.
In engineering practice, many advanced machine structures, such as the advanced aerospace crafts and shuttles, exposed to superhigh temperatures and the temperature or stress distribution could be in two or even three directions [12]. erefore, conventional FGM may not be useful in the design of such structures as conventional FGM has variation of material properties in only one direction.
e mechanical behavior of multidirectional functionally graded material structure is a valuable research topic. For instance, Lü et al. [13] adopted the state spacebased differential quadrature method to analyze the bending and thermal deformations of BDFG beams in which Young's modulus varies exponentially in the thickness and longitudinal directions. Şimşek [14] investigated the free and forced vibration of BDFG beams under moving load. Dynamic equations of beams in matrix form were obtained using Lagrange equations. Parametric studies of the fundamental frequency and dynamic response of BDFG beams were presented for different boundary conditions. Deng and Cheng [15] deduced the motion equations of the BDFG Timoshenko beam using Hamilton's principle. e dynamic stiffness method and modal superposition method were applied to calculate the free and forced vibration characteristics of beams, respectively. Fang et al. [16] researched the three-dimensional free vibration and time response of rotating BDFG cantilevered beams, where the material properties of functionally graded beams are assumed to change gradually through both the width and the thickness in power-law form. Based on the quasi-3D shear deformation theory, Karamanli [17] investigated the static bending behavior of BDFG sandwich beams with different boundary conditions under the action of uniformly distributed load. e displacement and stresses are solved using the Symmetric Smoothed Particle Hydrodynamics method. Wang et al. [18] presented the analytical solutions of the natural frequencies of a type of BDFG beams in which the material properties vary exponentially in axial direction and have a power law gradation in thickness direction. Truong et al. [19] presented an artificial neural network-differential evolution method to optimize the material distribution in BDFG beam under free vibration. In another work, Truong et al. [20] investigated the material optimization problem for BDFG beam under static load using the differential evolution method. Fariborz and Batra [21] studied the free vibration of curved BDFG circular beams using shear deformation theory and explored the mode veering phenomenon. Lei et al. [22] studied the postbuckling response of BDFG imperfect beam, where the critical buckling loads and the postbuckling equilibrium paths are under different boundary conditions using the generalized differential quadrature method. Tang et al. [23] studied the linear and nonlinear free vibration of BDFG beams with different boundary conditions. ey obtained the natural frequencies and mode shapes of BDFG beams using the generalized differential quadrature method and found that 2D FGMs may induce asymmetric modes in free vibration. e closed-form solutions of nonlinear frequencies of BDFG beams are determined by combining the Galerkin technique and the homotopy analytical method. In another work, Tang and Ding [24] analyzed the nonlinear free vibration of BDFG beams under hygrothermal effects. Recently, the small scale effect on the mechanical behavior of BDFG beams has drawn attention of researchers. For example, Chen et al. [25] analyzed the linear buckling, vibration, and stability behaviors of BDFG microbeams embedded in elastic medium using the modified couple stress theory. Considering the von Karman nonlinear strain, Chen et al. [26] investigated the nonlinear large deflection of BDFG microbam in the postbuckling domain, and the mode veering phenomenon is explored in the postbuckling vibration. Based on the modified couple stress theory accompanying with transverse shear-normal deformation beam theory, Karamanli and Aydogdu [27] investigated the free oscillation of the rotating BDFG porous sandwich microbeams and found that the lightweight and low-cost objectives can be achieved by controlling the material and porosity distribution through the structure. Yu et al. [28] studied the bending and free vibration of small-size BDFG beam by employing NURBSbased isogeometric analysis combined with a nonclassical quasi-3D beam theory. Bhattacharya and Das [29] investigated the nonlinear chord-and flapwise vibration of BDFG double-tapered rotating microbeam by considering the high operating temperature. Rajasekaran and Khaniki [30] obtained the forced vibration response of nonuniform BDFG microbeam under a moving harmonic load/mass. e dynamic reaction of microbeam is acquired by adopting the implicit finite element methods combined with various steps and Wilson-theta method. Adopting the nonlocal strain gradient theory, Sahmani and Safaei [31] applied different homogenization schemes to analyze the postbuckling and nonlinear bending behavior of BDFG micro/nanobeams. Considering the thermal effect, Karami et al. [32] deduced the tapered BDFG porous microbeam model based on nonlocal strain gradient theory. e free vibration problem is solved using the generalized differential quadrature method. Faroughi et al. [33] studied the wave propagation of BDFG porous rotating nanobeams on the basis of general nonlocal theory and Reddy's thirdorder shear deformation theory. Considering the temperature-dependent material and nonlocal elasticity theory, Lal and Dangi [34] studied the thermomechanical vibration of BDFG nonuniform nanobeam under linear and nonlinear temperature profiles.
From the literature review presented above, one may realize that the previous works on the dynamic behavior of BDFG beams are mostly limited to linear free vibration. To the author's best knowledge, there is no investigation related to the nonlinear parametric dynamics of bidimensionally FG beam in the open literature. e main aim of this paper is to present the parametric resonance of BDFG beam under the action of a time-dependent axial load for the first time. To formulate an accurate dynamic mathematical model of beam, the coupled nonlinear partial governing equations retaining the transverse and longitudinal displacements and inertia as well as the rotation and rotary inertias are established using Timoshenko beam theory and von Karman's nonlinear theory is considered. e parametric resonance behavior is acquired employing a two-step solving strategy. First, using the Galerkin technique, the continuous model of beam is discretized into the reduced order model which is a high-dimensional system of ordinary differential equations capturing almost all modal interactions and internal energy transfer between different modes. en, the reduced order model is solved for the period motions of beam using pseudoarclength continuation method. In the numerical results part, the influence of system parameters on the parametric dynamic behaviors of BDFG beams is discussed. Particularly, the dominant role of material distribution pattern in the parametric dynamics of BDFG beam is highlighted.

2
Shock and Vibration Figure 1 depicts the BDFG beam subjected to axial loads P ex in a rectangular coordinate system. e geometric parameters of beam are length L, width b, and thickness h. e beam is made of at least two different materials and the material properties gradient changes along both thickness and axial directions. e effective material property P(x, z), such as Young's modulus E, shear modulus μ, and mass density ρ of beam, can be calculated using the rule of mixture [35,36]:

Mathematical Formulation
where P i and V i are the material properties and volume fractions of the constituents. In this research, the volume fractions follow the power law of distribution throughout the axial and thickness direction of beam in the form of . Four different material distribution patterns are given as follows: TypeI: TypeII: TypeIII: TypeIV: in which m and n are the material gradient parameters.
On the basis of Timoshenko beam theory, the displacement field of BDFG beam can be written as where u 1 , u 2 , and u 3 are displacement components of any material point in x, y, and z directions, respectively. u and w are the axial and transverse displacements of the midplane along x and z directions. ψ is the rotation of transverse normal at midplane.
Using the displacement field (6) and the classical elasticity theory, one obtains the nonzero strain components of beam as where ,x or ,z represents the first-order derivative with respect to the coordinate x or z. Moreover, the constitutive equations are expressed as follows: Hamilton's principle [37] is employed to derive the governing equations of BDFG beam, which reads e first variation of strain energy is expressed as in which the forces and moments resultants are defined as where k s � 5/6 is the shear correction factor. Moreover, the stiffness terms are given as e first variation of kinetic energy of beam is where the mass inertia terms are defined as Shock and Vibration e virtual work done by the axial force P ex and the damping force can be written as where c d and c r are the damping coefficients. Moreover, the axial force P ex is in the form of in which P s and P d are the static and dynamic amplitudes of axial load.
Substituting (10), (13), and (15) into (9) and introducing the dimensionless quantities the nondimensional governing equations can be obtained as in which the dimensionless coefficients are given as Furthermore, the dimensionless simply supported condition is given as

Solution Method
e nonlinear partial differential equations with variable coefficients presented in (18)-(20) are continuous system with infinite degrees of freedom. Generally, such a continuous system is intricate to solve directly. Krysko et al. [38,39] employed the finite difference method and finite element method to solve the nonlinear mathematical models of beam (without variable coefficient) for the time histories of response. However, this strategy is heavily dependent on the numerical value precision and not applied well to bifurcation analysis of periodic response. A different tack is discretizing the continuous model to lowdimensional system and then deals with approximate analytical methods [40][41][42][43][44][45]. Practically, the features of nonlinear vibrations and dynamics of a continuous system are described by too many degrees of freedom. In this study, each equation of motions given in (18)-(20) is transformed into high-dimensional nonlinear ordinary differential equation system, namely, the reduced model. e following series expansions are considered as the approximate solution of motions: where p m (t), q m (t), and r m (t) are the generalized coordinates of the longitudinal, transverse, and rotational motions, φ m (x) and ϕ m (x) are the eigenfunction for the longitudinal and transverse vibration of a uniform simply supported beam, and ξ m � ϕ m, x /(mπ). Inserting expressions (23) into (18)- (20) and using Galerkin's method, one obtains In the practical calculation, M � N � Q � 8 is chosen such that the results are converged. e linear eigenvalue problem is easily solved by dropping the damping and nonlinear terms in the previous system. For nonlinear analysis, the response of system is obtained by direct time integration. e periodic solutions of equations (24)- (26) are traced numerically utilizing pseudoarclength continuation technique [46,47]. Specifically, the frequency-and force-response curves of beam are achieved by using a numerical continuation procedure.

Numerical Examples
In this section, numerical results for primary resonant response of simply supported BDFG beam are provided. It is assumed that the BDFG beam is made of a mixture of ceramic (SiC: E 1 � 427 GPa, ρ 1 � 3100 kg/m 3 , and ] 1 � 0.17) and metal (Al: E 2 � 70 GPa, ρ 2 � 2702 kg/m 3 , and ] 2 � 0.3). e width and the thickness of the BDFG beam are taken to be equal to each other (b � h). In the following analysis, numerical results are evaluated at the specified point (x � 0.6) of the beam.
To check the accuracy of the present mathematical model, we present a comparison of the dimensionles fundamental frequency of a simply supported BDFG beam between the present results and those of Şimşek [14], as shown in Table 1. e material properties and geometry of the BDFG beam given in [14] are used. Table 1 demonstrates that the present results are very close to those given by Şimşek [14]. Further, it demonstrates that the present mathematical modeling is suitable for different type of material distribution in two dimensions.
In Figure 2, the effects of gradient parameters on the dimensionless fundamental frequencies of BDFG beams are examined by considering different type of material distribution, when P s � 5 and P d � 0. It is notable that the dimensionless fundamental frequencies will descend with the increase of gradient parameters m and n. is descent is due to the fact that the volume fractions of ceramics V c will diminish as m and n get larger (as shown in equations Shock and Vibration (2)-(5)), which leads to the reduction of the elasticity modulus. Another important observation is that the downtrend of dimensionless fundamental frequency caused by the gradient parameters will be more remarkable, when m is small than 5. Further, the type of material distribution of BDFG has important influence on the frequency, and the order is III > I > IV > II.
In the following part, we explore the principal parametric resonance behaviors of BDFG beams by performing the bifurcation analysis of the periodic motion of system. e frequency-force-response curves and time histories are plotted to examine the effects of system parameters on the nonlinear parametric dynamics of beams. Unlike the primary resonance of BDFG beams [48] which is transversely excited and the resonance response which appears, when being near the linear fundamental frequency ω, the principal parametric resonance of an axially excited BDFG beam occurs, when the excitation frequency Ω is nearly double linear frequency (2ω l ). In Figure 3, the frequency-response curves for the transverse, longitudinal, and rotational motions of type I BDFG beam are plotted, when m � n � 1, L/h � 20, and c d � c r � 0.05. In this case, the static axial load is P s � 5 and the amplitude of dynamic load is P d � 1.5. For the chosen system parameters, the linear frequency of BDFG beam is ω 1 � 12.0059. As can be seen from Figure 3, the parametric dynamics have a trivial solution throughout the state space. As excitation frequency Ω approaches 2ω 1 , the branch point (BP) of cycle occurs at Ω � 23.6816 (BP 1 ) and 24.2993 (BP 2 ), where the trivial solution loses stability and nontrivial solutions appear. e results demonstrate that all the transverse, longitudinal, and rotational motions of beams exhibit hardening-type nonlinearities and cyclic-fold (CF) bifurcation points are detected at Ω � 27.8025, where the jump phenomenon happens. Figure 4 depicts the frequency-response curves for the transverse motions of BDFG beam with different type of material distribution, when m � n � 1. e geometric parameters and the external load used in this figure are the same as in Figure 3. It is noted that the type of material distribution has great influence on the parametric resonance dynamics. More specifically, the variation of material distribution mode will change the width and location of resonance region, the maximum amplitude of period motion in resonance region, and the breadth of unstable trivial solutions. In Figure 5, the time histories and phase diagrams of the transverse motions of BDFG beams are displayed, when Ω � 2ω 1 . From Figures 4 and 5, a significant observation form is that the responses of type I and III BDFG beam are asymmetric about the undeformed static equilibrium position, but the responses of type II and IV beam are symmetric.
is unsymmetrical appearance is due to the stretching-bending coupling that results from the asymmetry of the material distribution in thickness direction. In type I and III BDFG beam with stretching-bending coupling (i.e., B xx ≠ 0), the strain energy contains odd powers terms such as ψ ,x w 2 ,x , which makes that the absolute values of maximal and minimum deflections on the two sides of static equilibrium position are not equal, when the strain energy is constant. Figures 6 and 7 examine the effect of axial load including static and dynamic parts (P s and P d ) on the frequencyresponse of transverse motion of type I BDFG beam. As shown in those two figures, the existence scopes of nontrivial solutions and the ranges of unstable trivial solutions are widened as P s and P d increase. Moreover, the raise of P s and P d results in the CF point occurring at large value of Ω. By increasing P s , the unstable trivial solutions shift left, but the change of P d will not affect the position of the region of unstable trivial solutions.
is is due to the fact that increasing the static part of axial load, P s , will lead to the softening of beam, whereas dynamic part P d is a periodic force which will not change the equivalent stiffness of beam.
To illustrate the influences of gradient parameters on the parametric dynamic of beam, we fix one gradient parameter (m or n) and let another vary; the results for the frequencyresponse curves of the type I BDFG beam are plotted in Figure 8. As can be seen from this figure, the rise of functional gradient parameters leads to a softer structure. More specifically, the parametric resonance region gets wider and shifts to left and the maximized amplitude of the period motion in resonance region gets larger substantially, when m and n increase. e reason for this phenomenon is that the increase of gradient parameters m and n will reduce    , one may realize that the effect of m (in axial direction) on the frequency-response curve is much significant than that of n.
Next, by fixing the excitation frequency Ω at a special value (e.g., the vicinity of twice the linear fundamental frequency) and choosing the magnitude of dynamic part P d as the variable parameter, the parametric force-response curve of BDFG beam is plotted in Figure 9. In this figure, the frequency of P d is fixed at Ω � 1.01 × 2ω l . e results show that the trivial solution loses stability via a BP point at P d � 1.4456, where the nontrivial solution occurs. e occurrence of BP point is fairly critical from the viewpoint of engineering, due to the fact that the response of system jumps from zero to a relatively large amplitude at this point as P d increases. A CF point is detected at P d � 1.32126 at which the nontrivial solution loses stability and the amplitude of periodic response drops swiftly to zero. Figure 10 depicts the force-response curves for different value of excitation frequency, e.g., Ω/(2ω l ) � 1, 1.01, 1.02, 1.03. As shown in this figure, the BP point shifts to right as Ω increases. However, the variation of Ω changes the position of CF point slightly. Particularly, there is no CF bifurcation, when Ω � 2ω l .
In Figure 11, the differences of the force-response curves of BDFG beam with various type of material distributions are explored, when Ω � 1.02 × 2ω l . It is remarkable that the type of material distribution has a significant effect on the force-response curves of BDFG beam. For example, the values relevant to bifurcation points (including BP and CF points) of type III beam are much larger than those of type I beam. Also, the amplitude of the response of type III beam is larger than other types of BDFG beam. Taking the type I beam, Figure 12 highlights the influence of gradient parameters on the force-response curves.
e results demonstrate that the bifurcation points occur as smaller value of P d as m and n get larger. Moreover, the effect of gradient parameter m on the position of bifurcation points is more remarkable than gradient parameter n; in contrast, n has more significant effect on the amplitude of response than m. e effect of the damping coefficients on the parametric      dynamic of type I BDFG beam is elucidated in Figure 13. From Figure 13, one can notice that the larger values of c d lead to the response amplitude decreases and the range of unstable trivial solution is wider. For the force-response curves shown in Figure 13(b), the increase of damping coefficients will make the BP and CF points shift left and the jump phenomenon will occur earlier.

Conclusions
e principle parametric resonance of bidirectional functionally graded beam is investigated. e BDFG beam is subject to a variable axial load. e material properties of beam vary in thickness and axial directions and four types of material distribution are considered. e continuous model of beam is formulated by means of Hamilton's principle in conjunction with Timoshenko beam theory.
e continuous system is discretized into a reduced order model employing the Galerkin technique. Principle parametric resonance behaviors of beam are simulated by preforming bifurcation analysis on the reduced system by the method of pseudoarclength continuation technique. e numeric results are depicted to illustrate the influence of system parameters. e main results of this paper are included as follows. All the transverse, longitudinal, and rotational motions of BDFG beam exhibit hardening-type behavior. e cyclic-fold bifurcations and branch point are detected in the frequency/force-response curves, and those two types of bifurcation may lead to the occurrence of jump phenomenon. e variation of material distribution will not change the type of nonlinearity and bifurcations but give rise to the quantitative difference of results. e asymmetrical material distribution in thickness direction of beam leads to the asymmetry of dynamic responses.

Data Availability
e data used to support the findings of the study are available upon request from the corresponding author.

Conflicts of Interest
e authors declare that there are no conflicts of interest regarding the publication of this paper.