Nonlinear Forced Vibration of Bidirectional Functionally Graded Porous Material Beam

The nonlinear forced vibration of bidirectional functionally graded porous material beams where the material components gradient change in both thickness and axial directions are studied in this study. Combining von Karman’s geometric nonlinearity and ﬁrst-order shear deformation theory, the governing equations describing the coupled deformations are formulated as a system of nonlinear partial diﬀerential equations. Utilizing the Galerkin method, the formulated continuous model is transformed to a coupled nonlinear ordinary diﬀerential dynamic system. By accomplishing bifurcation calculation for periodic response of the discrete system using pseudoarclength technique, the vibration response curves are obtained by extracting the max-min amplitude of periodic motions. To highlight the eﬀect of nonlinearity, the linear and nonlinear dynamic responses of beam are demonstrated. It is found that the periodic motion of beam may undergo cyclic-fold bifurcation. Numerical results are presented to examine the eﬀects of the system parameters, e.g., gradient indexes, porosity, damping coeﬃcients, and aspect ratio.


Introduction
Functionally graded materials (FGMs) are a class of composites which have a continuous variation of material properties in one or more directions and hence have the ability to lessen the phenomenon of stress concentration that is usually unavoidable in laminated composites. Besides, FGMs can be designed according to the service condition, such as extreme high temperature or corrosive environment. With these merits, FGMs has been widely used in many industrial fields such as aircraft industry, aerospace manufacturing, and the atomic engineering. In the past few decades, static and dynamic behaviors of classical 1D-FGM beams, in which the grading of materials only exist in transversely or axially direction of beam, were studied extensively. For example, Kadoli et al. [1] investigated the bending response of metal-ceramic FG beams in thermal environment using higher order shear deformation beam theory. Niknam et al. [2] provided the analytical solution of the nonlinear bending of tapered functionally graded beam under action of thermal and mechanical load. Adopting the differential quadrature method (DQM), Esfahani et al. [3] performed thermal buckling analysis of FG beams sited on the elastic foundation with nonlinear reaction force. Considering the materials grading in thickness or axial direction, Alshorbagy et al. [4] showed the vibration characteristics of FG beams employing the finite element method (FEM). Şimşek and Kocatürk [5] analyzed the natural frequency and forced vibration response of FG beam under the action of a concentrated moving harmonic load. Shahba et al. [6] examined the frequency characteristic and dynamic stability of tapered FG beam using FEM. Atmane et al. [7] analyzed the free vibration of an exponential tapered FG beam, and the closed-form solution of frequencies is obtained. Some researchers are also interested in the nonlinear behaviors of functionally graded beam. For example, Ke et al. [8] gained the approximate analytic solution of nonlinear vibration frequency of FG beams with different end supports. Employing He's variational method, the nonlinear free vibration analysis of FG beams in thermal environment was presented by Fallah and Aghdam [9,10]. Shooshtari and Rafiee [11] adopted the multiple time scale method to study the forced nonlinear primary, superharmonic, and subharmonic resonances of FGM beams with immovable ends.
In modern industry, many advanced structures are exposed to complicated service environment which leads the temperature or stress distribution grading in several directions [12]. is reality makes the traditional FGMs not applicable in those structures on account of material component that only vary in single direction. us, the multidirectional FGMs structure is a worthy research theme. As example, Fang et al. [13] analyzed the vibration frequency of rotating BDFG cantilevered beam. In their beam model, the materials properties of beam change gradually in the width and thickness. Karamanlı [14] formulated an investigated quasi-3D shear deformation beam model to study the bending of BDFG sandwich beams. Truong and co-workers [15,16] proposed an evolution method for material optimization in the design of bending and free vibration behaviors of BDFG beam. Fariborz and Batra [17] researched free vibration of curved BDFG circular beams and explored the mode veering phenomenon. Using the DQM, Lei et al. [18] analyzed the critical buckling load and static response of BDFG imperfect beam subjected to axial load. Tang and Ding [19] investigated the nonlinear vibration frequency of the BDFG beams considering the hygrothermal effects.
e FG porous beam also attracts researchers' attention. Chen et al. [20] investigated the linear vibration response of FG porous beams. Results show that the symmetric porosity distribution of porosity gives the highest vibration frequency and, meanwhile, the smallest response amplitude. In another work, Chen et al. [21] investigated the large amplitude free vibration of sandwich porous FG beam. Ebrahimi and Zia [22] acquired the approximate expression of nonlinear vibration frequency of FG porous beam applying the multiple scale method. On the basis of fully geometrically nonlinear beam theory, Tian et al. [23] formulated the large deformation vibration of rotating double-tapered FG porous beam considering the interaction among bending, stretching, and twist modes. Employing the finite element method and plane solid continua model, Akbaş [24] analyzed the linear vibration response of FG porous deep beam and showed that the porosity has obvious significant influence on the dynamic responses of beam.
From the above literature review, one may draw a conclusion that open literatures on the kinetic characteristics of BDFG beams are mainly focusing on the linear vibration behaviors. In the open literature, the nonlinear forced vibration response of BDFG porous beam with different boundary conditions has not been reported. In this research, we attempt to explore the nonlinear forced vibration response of BDFG porous beam under action of lateral harmonic excitation. e mathematical model is established based on the Timoshenko beam theory in conjunction with von Karman's nonlinear strain.
en, the reduced order model is obtained by employing the Galerkin procedure. e periodic motion of BDFG porous beam is acquired by dealing the reduced order model with the pseudoarclength method. Numerical simulations are given to explain the influences of material, porosity, and geometrical parameters on the nonlinear primary resonance of BDFG beam. Figure 1 depicts the BDFG porous beam with following geometric dimensioning. e constitutes of beam are twophase materials. e material properties' gradient changes in axial and thickness directions of beam. e effective material properties of beam can be estimated employing the mixture rule [25]:

Modeling
where P i (i � 1, 2) is the material properties of the constituent, and α(x, z) is defined as the porosity volume fraction ratio. V i is the volume fraction of the constituent and satisfy the following relation [26,27]: where m x and m z represent the gradient varying indexes. Moreover, α(x, z) is defined as where α 0 is the porosity volume fraction coefficient. e displacement field of Timoshenko beam is given as where u i with i � 1, 2, 3 is the displacement component of an arbitrary material point of beam in space coordinates. u, w, and ψ are the stretch, transverse, and shear deformation components of a point at the mid-plane of beam. Employing the described displacements above, the strain components of Timoshenko beam are where ,x and ,z represent the first-order differential with respect to space coordinates x and z. e constitutive relations are given as where μ(x, z) � (E(x, z)/2(1 + ](x, z))). e Hamilton's principle [28] will be applied to formulate the dynamical model of BDFG beam. at is, e first variation of strain energy of beam can be represent as 2 Shock and Vibration in which, where k s is the shear correction factor in Timoshenko beam theory. Besides, the rigidity coefficients are calculate as e variation of kinetic energy is given as in which, the inertia coefficients are e virtual works introduced by the external excitation and damping are in which, Ω is the frequency of external excitation; c d and c r are the damping coefficients. Substituting equations (10), (11), and (13) into equation (9), one has the following equilibrium relation in the form of forces and moments: Defining the following nondimensional variables, e dimensionless nonlinear motion equations are gained as where the nondimensional coefficients are defined as In addition, the dimensionless boundary conditions for simply supported (S) and clamped (C) ends are given as

Solving Method
In this study, the approximate solutions of the continuous system in equations (16)- (18) are assumed as where p i , q i , and r i denote the modal coordinates of the coupled motions including the longitudinal, transverse, and rotational deformations, respectively. Moreover, φ i is the mode function of longitudinal vibration and ϕ i is the eigenfunction of transverse vibration of uniform beam. ξ m is the first-order differential of mode function of transverse vibration. e external excitation is assumed as F � fϕ 1 (x).
Substituting the approximate solution equation (21) into equations (16)-(18) and employing Galerkin's method, we have the reduced order model as follows: It is easy to perform linear vibration analysis by eliminating the damping, nonlinear, and excitation terms in the above discrete system. e nonlinear vibration response of beam can be acquired by solving equations (22)-(24) with the time integral method. e periodic solutions of nonlinear vibration are traced numerically by the means of pseudoarclength continuation technique [29], and then, the maximum and minimum values of periodic motions are extracted to construct the frequency-force response curves. Concretely, the response curves of beams are achieved with two steps: (i) first, the frequency of external excitation Ω is fixed far away from the linear vibration frequency of BDFG beam, and f is set as the bifurcation parameter. e solution is started from zero solution, and then, f is increased to the desired value. (ii) As f reaches the desired value, the solution is continued by choosing Ω as the bifurcation parameter.

Results
In this part, numerical results and discussion for nonlinear vibration of BDFG porous beams are presented to analyze the influence of system parameters. e various end supports including the double clamped (CC), clamped-simply supported (CS), and double simply supported (SS) edges are considered. e material components of beam are ceramic (SiC: E 1 � 427 GPa, ρ 1 � 3100 kg/m 3 , ] 1 � 0.17) and metal (Al: E 2 � 70 GPa, ρ 2 � 2702 kg/m 3 , ] 2 � 0.3). In the following results analysis, periodic responses of specified point (x � 0.6) are used to evaluate the response curves of beam.
To verify the derived mathematical model, the present dimensionless fundamental frequencies are contrasted to those given by Şimşek et al. [30], as tabulated in Table 1. For comparison, the material distribution pattern is defined as P(x, z) � P 0 e k 1 (x/L+0.5)+k 3 (z/h+0.5) , where P 0 is the reference material parameter, including E 0 � 210 GPa, ρ 0 � 7850 kg/m 3 , and ] � 0.3. e results illustrated in Table 1 demonstrate that a good agreement is achieved. For further validation of the present nonlinear vibration response, we presented a comparison of the frequency-response of simply supported BDFG beam in Figure 2. As interpreted in this figure, the presented frequency-response curve is perfectly in tune with those given in Ref. [31]. In Figure 3, the convergence analysis of the reduced order model is executed when m x � m z � 1, f � 0.2, and c d � c r � 0.05. As shown in this figure, the results will be obtained reasonably when the mode number is equal to 6. In the following calculation, P � Q � R � 8 is used to guarantee the convergence of results.
In Figure 4, frequency-response curves of the SS BDFG porous beam are presented when m x � m z � 1, α 0 � 0.1, L/h � 20, c d � c r � 0.05, and f � 0.2. is figure shows that the nonlinear vibration of beam demonstrates hardening-type nonlinearity. Varying the excitation frequency Ω, periodic motion of beam will give rise to cyclic-fold (CF) bifurcation, and two CF points are spotted in the presented primary resonance region. At the bifurcation point, the nonlinear vibration response of the system leaps from one stable periodic solution to another. Moreover, one can recognize from Figures 4(a) and 4(b) that the maximum and minimum values of the periodic solution are unsymmetrical with respect to trivial undisturbed configuration of beam. e differences of linear and nonlinear dynamics of SS BDFG porous beam are explored in Figure 5. In Figures 5(a), 5(d), and 5(g), the linear and nonlinear response curves depicting the lateral motion of BDFG beam are demonstrated when m x � m z � 1, α 0 � 0.1, L/h � 20, c d � c r � 0.05, and f � 0.2. It can be noticed that the maximum amplitude of linear response is larger than that of the nonlinear response. Comparing to the linear BDFG porous beam in which the resonant region is located nearby the linear frequency, resonant region of nonlinear BDFG porous beam shifts right due to the hardening-type nonlinearity. e time histories (Figures 5(b), 5(e), and 5(h)) and phase diagrams ( Figure 5(e), 5(f ), and 5(i)) of linear and nonlinear responses are illustrated when Ω � 14. It is observed from these subgraphs that all the linear and nonlinear responses of transverse and rotational deformations exhibit smooth periodic motions with similar features. However, the linear and nonlinear responses of longitudinal motions display considerably difference. Specifically, the smooth periodic solution is presented in linear longitudinal motions. For nonlinear response, conversely, a small snap through motion is observed.
In Figure 6, the frequency-response diagrams for different value of m x and m z are plotted to explore the effects of functional gradient indexes on transverse motion of SS BDFG porous beam. As shown in Figure 6, increasing the functional gradient parameters may result in the softer BDFG beam. More specifically, the resonance region turns to smaller value of Ω and the magnitude of periodic motion increases substantially as m x and m z get larger. From Figures 6(a) and 6(b), one can recognize that the effect of gradient parameter m x (in axial direction) on nonlinear vibration features of BDFG porous beam is much significant than that of m z .
In Figure 7, the force-response curves of SS BDFG porous beam are depicted when the excitation frequency Ω � 15. As shown in this figure, two CF points are spotted. As f varies, the periodic solution of BDFG porous beam undergoes a leap at points CF 1 resulting in other periodic attractors with larger amplitude. When f decreases from 1, another jump occurs at points CF 2 resulting in smaller periodic motion. Here again, the dynamic responses of the system are asymmetrical (Figures 7(a) and 7(b)), caused by asymmetric distribution of the material. e force-response curves of SS BDFG porous beam are displayed for various excitation frequencies as shown in Figure 8. e geometric and material parameters are same to those in Figure 7 except excitation frequency. As shown in this figure, CF points occur for some specific excitation frequencies (e.g., Ω � 14.5, 15, and 15.5). Furthermore, the coexistence region where two stable periodic solutions coexist between two CF points will shift to higher value of f and get broader as the increase of Ω. Figure 9 inspects the effects of graded parameters on the force-response curves of SS BDFG porous beam. e excitation frequency Ω � 1.1ω l . Here again, increasing m z (or m z ) gives rise to the amplitude of motion. However, the larger values of m z and m z result in smaller coexistence region in which two stable periodic solutions coexist. Shock and Vibration Figure 10 depicts the variation of frequency-response curves as the magnitude of excitation f and damping coefficients change when m z � m x � 1 and L/h � 20. It is observed from Figure 10(a) that the larger values of damping coefficients lead to the narrower primary resonance region. e maximum amplitude of the response will increase as c d and c r decrease. Instead, as shown in Figure 10(b), a larger f leads to wider resonance region and larger response amplitude. e frequency-response curves of the nonlinear transverse motion obtained for different L/h are illustrated in Figure 11.
e system parameters are chosen as m z � m x � 1, f � 0.2, and c d � c r � 0.05. In numerical     e results show that the resonance region gets wider, but the maximum amplitude of the response in the resonance region will decrease when aspect ratio L/h increases.
us, one may draw the conclusion that the hardening-type nonlinearity is much prominent for a slender BDFG porous beam. Figures 12 and 13 demonstrate the influence of porosity including the volume fraction coefficient α 0 and the distribution pattern of porosity on nonlinear vibration response of BDFG beam. As shown in Figure 12, the increase of α 0 will shift the resonance region of beam to lower frequency in both types of porous distribution. e maximum amplitude of response in the resonance region of nonuniform porosity beam will minish as α 0 increases, but this changing trend is inverse in uniform porosity beam. Besides, the influence of α 0 on the frequency-response curve of uniform porosity beam is much significant than that on nonuniform porosity beam. e same phenomenon is also observed in Figure 13. From this figure, one may deduce that the larger α 0 leads higher amplitude of the force-response curve. To probe the effect of boundary conditions on vibration response, Figure 14 plots the frequency-and force-response of BDFG porous beams with SS, CS, and CC edge supports. As shown in Figure 14(a), the resonant frequency region of CC BDFG porous beam occurs at the highest value of Ω and that of SS BDFG beam is the lowest. It should be pointed out that the maximum amplitude of the response in the resonant region of CS BDFG porous beam is higher than that in SS BDFG porous beam.
is phenomenon is owing to the fact that the chosen sampling point, x � 0.6, is close to the maximum lateral displacement position of CS BDFG porous beam and deviate   from that of SS BDFG porous beam. From Figure 14(b), one can realize that the relation of the value of f corresponding to the occurrence of CF points is CC > CS > SS. Moreover, the width of coexistence region obeys the same pattern.

Conclusions
is research focuses on the nonlinear forced vibration of BDFG porous beam. e material composition changes gradually in thickness and axial directions of beam complying with a power-law. e motion equations of beam are formulated by using Timoshenko beam theory. e reduced order model is obtained by discretizing the original system using the Galerkin procedure. Nonlinear vibration responses of the reduced system are acquired by the means of numerical integration and pseudoarclength technique. e primary outcomes of this study are listed as follows: (i) the couple nonlinear responses of BDFG porous beam reveals hardening-type behavior, and the periodic response of beam undergoes cyclic-fold bifurcations as frequency and amplitude of external excitation vary; (ii) for considering the material distribute pattern, the BDFG porous beam gets softer with the increase of functional gradient parameters; (iii) the resonant region shifts to lower excitation frequency and gets broader as functional gradient parameters increase; (iv) increasing the amplitude of the excitation (or decreasing the value of the damping) results in a wider response region and larger amplitude of response; (v) the porosity distribute pattern shows an appreciable impact on vibration response of beam; and (vi) the hardening-type nonlinear behaviors of BDFG porous beam become more remarkable as aspect ratio increases.

Data Availability
No data were used to support this study.

Conflicts of Interest
e authors declare that they have no conflicts of interest.