Abstract
This paper presents a detailed study on the stochastic stability, jump, and bifurcation of the motion of the composite laminated beams subject to axial load. The largest Lyapunov exponent which determines the almost sure stability of the trivial solution is quantificationally resolved and the results show that the increase of the bandwidth facilitates the almost sure stability of the trivial response. The stochastic jump and bifurcation of the response are numerically calculated through the stationary joint probability and the results reveal that (a) the higher the excitation frequency is, the more probable the jump from the stable stationary nontrivial solution to the stable stationary trivial one is; (b) the most probable motion is around the nontrivial solution when the bandwidth is smaller; (c) the outer flabellate peak decreases, while the central volcano peak increases as the value of the excitation load decreases; and (d) the overall tendency of the response is that the probable motion jumps from the stable stationary nontrivial branch to the stable stationary trivial one as the fiber orientation angle of the first lamina with respect to the -axis of the beam increases from zero to a smaller angle.
1. Introduction
Laminated composite beams are basic structural elements and widely used in many engineering fields such as large space station, aircraft, automobiles, high-speed mechanism, high-pressure vessel, and submarine. In fact, these laminated beams often operate in complex environmental conditions and are exposed to abominable dynamic excitations resulting in excessive vibration and fatigue damage. In the case of moderately thick laminated beams, the classical laminated beam theory can become inaccurate due to neglecting the transverse shear and normal strains in the laminate [1]. In order to take the effect of low ratio of transverse shear modulus to the in-plane modulus into consideration, the first-order shear deformation theory, which gives emphasis on the constant assumption of the transverse shear strain across the thickness of the beam, has been developed [2, 3]. Compared to the first-order shear deformation theory, the higher-order shear deformation theory does not require a shear correction coefficient and gives more realistic representation of the transverse deformation vibration. A number of higher-order theories with different shear strain shape functions (including polynomial functions [4–8], trigonometric functions [9–11], and exponential functions [12, 13]) have been proposed.
Naturally, some researchers focus on the dynamics of composite laminated beams based on various theories. Kant et al. [7] proposed an analytical method for the dynamic analysis of laminated beams using higher-order refined shear deformation theory. Song and Waas [14] investigated the bucking and free vibration of stepped laminated beams using a simple higher-order shear deformation theory. Matsunaga [8] studied the natural frequencies, buckling stresses, and interlaminar stresses of composite beams with simply supported edges using a global higher-order theory. Subramanian [15] performed the free vibration analysis of laminated composite beams using two higher-order displacements, that is, a quintic and a quartic variation of in-plane and transverse displacements in the thickness coordinates of the beam.
Out of question, the nonlinear dynamic research on the composite laminated beams is far from abundance. Over the past few decades, great attention has been paid to the analysis of the complex nonlinear dynamics of isotropic, flexible beams. Crespo Da Silva and Glynn [16, 17] found that the generally neglected nonlinear terms due to curvature are of the same order as the nonlinear terms due to inertia and that the curvature terms may have a significant influence on the response of the inextensional flexural-flexural-torsional beams. Nayfeh and Pai [18] and Pai and Nayfeh [19] used the equation of motion formulated in [16] to analyze the nonlinear vibration of a cantilever beam subject to principal parametric and primary excitations and found that the geometric nonlinear terms have a hardening effect, whereas the inertia terms have a softening effect. Zavodney and Nayfeh [20] derived the nonlinear partial differential equation for a slender cantilever beam carrying a lumped mass at an arbitrary position and investigated the principle parametric resonance of the single mode in both theory and experiment. Following the equation of motion of the beam described in [20], Kar and Dwivedy [21] and Dwivedy and Kar [22–26] systematically dealt with the nonlinear dynamic behaviors of a slender beam carrying a lumped mass with principal parametric, combination parametric, and internal resonance of the lower modes. Anderson et al. [27] experimentally investigated the nonlinear resonances of a flexible beam among its first four natural frequencies subject to external or/and internal excitation. Anderson et al. [28] focused on the investigation of the nonlinear coupling behaviors of a thin, slightly curved, isotropic, flexible cantilever beam between its high-frequency modes and a low-frequency mode in both theory and experiment; they developed an analytical model to explain the interactions between the widely separated modes and used a three-mode Galerkin truncation to obtain a sixth-order nonautonomous system. Particularly, their analytical prediction on the responses and bifurcation diagrams are in good qualitative agreement with the experimental observations.
In most studies of the nonlinear dynamics of aforementioned beams, the excitation is a pure harmonic one. Naturally, some problems still need to be elaborated. For instance, Zavodney and Nayfeh [20] used a signal generator and a power amplifier to drive a modal shaker to produce base excitation so as to capture experimentally the frequency-response curves of the first modal principal parametric resonance. Results in Figures 10 and 16 of [20] show that the frequency-response curve is a definite overhang; that is, jumps from the higher branch (nontrivial) to the lower one (trivial) occur or vice versa, which does not accord with theoretical analysis. Anderson et al. [29] experimentally investigated the only planar response of a parametrically excited slender inextensional cantilever beam based on the analytical model in [16] using the similar excitation equipment as those in [20]. Their experimental results showed that the frequency-response curves for the first two modal principal parametric resonances are also in the form of definite overhang. They added a quadratic damping into the nonlinear dynamic equation and made the experimental and theoretical results be in agreement on both the frequency-response and force-response curves of the first two modes. Feng et al. [30–32] thought that it is very difficult to accurately keep the excitation frequency of the shaker at a fixed one to experience a long-term excitation. In other words, the excitation frequency will slightly drift off the center (fixed) frequency in random. Thus, they introduced the so-called narrow-band bounded noise to feature the shaker excitation and their research results show that the jumping phenomena in the experimental results of [20, 29] can be explained and described using the theory of random vibration. In other words, such experimentally investigated jumps in [20, 29] may be the stochastic jump in nature.
Since references to the nonlinear dynamics of composite laminated structures subject to random excitation are few up to now, in the present study, the stochastic principal parametric resonance of composite laminated beams is systematically dealt with. The nonlinear governing partial differential equations of the motion are derived by using higher-order shear deformation theory and Hamilton’s principle, which are suitable for composite laminated beam subject to the axial excitations. In particular, the excitation is modeled as a narrow-band bounded noise. Numerical calculations for the almost sure stability of the trivial solution and the stochastic jump and bifurcation of the response for the nontrivial solution are presented and results are discussed. The aim of the present work is to find the inherent characters of the composite laminated beams subject to narrow-band random excitation.
2. Formulation and Analyses
2.1. Equation of Motion
A simply supported laminated composite beam is shown in Figure 1 with length and thickness subject to axial load . The displacement proposed by Levinson [4] and used by Reddy [6] is given as belowwhere , and are total displacements, and are the displacements of a point in the middle surface of the beam in the and directions, and is rotation of the normal about the -axis.
According to Reddy and von Karman’s theory, the governing equations of the motion by the Hamilton principle can be obtained aswhere is the coefficient of damping and other coefficients can be found in the appendix.
Also, the boundary conditions can be written as : : and : and : and :
In what follows, we focus on the case of regular, symmetric cross-ply laminated composite beam. Thus, according to the character of the odd function, we have
Substituting (7) into (2a) and (2c) and neglecting the longitudinal inertia and the moment of inertia about the neutral axis, we have
Integrating (8) and using (3a) and (3b), we arrive at
According to the boundary conditions, the solution of (9) can be approximately expressed as [33]
In order to obtain the dimensionless equations, the variables and parameters are introduced as given below
Naturally, for simplicity, the star is dropped in the following analysis.
Closed form solutions for can be obtained for the beam by assuming
In what follows, we focus on the analysis of the principal parametric resonance of the first mode; it is also assumed that there is not any internal resonance between two modes. Therefore, out of the infinite modes present in and in the present of viscous damping, only the excited first mode will contribute to the long-term response. Thus, (2b) can be discretized by taking the first mode of (13) into consideration and using Galerkin’s method [33] and the nonlinear dynamic equation is finally obtained as where
In most studies of the nonlinear dynamics of aforementioned composite laminated beams, the excitation is a pure harmonic one; however, from the point of technical or engineering view, it is hard for a laminated structure to experience a pure and long-term harmonic excitation; in other words, such a harmonic excitation assumption is at variance with the reality. Following [30–32], in what follows, we use the narrow-band bounded noise to feature a kind of narrow-band random excitation; that is, where and are the amplitude and center frequency of the excitation, which are constants, is an intensity, which represents the bandwidth, of the narrow-band random excitation, is a standard Wiener process, and is a uniformly distributed random number in . The inclusion of the phase angle makes a stationary process. So is a non-Gaussian distributed stationary stochastic process and has density function , zero mean value, and spectral density function as follows:
The bandwidth of process depends mainly on parameter . It is a narrow-band process when is small and a wide-band process when is large. Thus, (14) can be rewritten as where , , , and is a small parameter.
Also, a new time scale is defined as . Let , , , and and rewrite (18) as
2.2. Analysis of Principal Parametric Resonance
In what follows, calculations have been performed for a graphite/epoxy symmetric laminated cross-ply beam with dimensions m, m, m, kg/m3, GPa, GPa, GPa, GPa, , and . Thus, parameter in (19) is calculated to be 2.4566 when the fiber orientation angle of the first lamina with respect to the -axis of the beam is zero.
In order to analyze the solution of the nonlinear equation subject to principal parametric resonance, the method of multiple scales is employed here. Thus, a first-order uniform expansion of the form is used as given below
Moreover, for unit Wiener progress , because of and , we have
Here the discussion and investigation are restricted to the case of principal parametric resonance of the first mode. Thus, the frequency detuning parameter is introduced to describe the closeness to the principal parametric resonance as given below Substituting (20) and (22) into (19) yields order : order :
The first-order approximate solution of (23) is where cc stands for the complex conjugate of the preceding terms.
Substituting (22) and (25) into (24) and eliminating the secular terms, we have where .
Substituting the polar form into (26) and separating the real and imaginary parts, we get
2.2.1. Stability Analysis of the Trivial Response
It can be found that (28a) and (28b) have a trivial steady state solution . In order to determine the stability of the trivial steady state response, the trivial response is expanded near and the corresponding nonlinear terms are neglected; the linearization form of (28a) and (28b) can finally be expressed as
Let and we rewrite (29a) and (29b) in the form of Itô equation as given below
It can be found that the stochastic process generated on by (30a) and (30b) is Markov and is also ergodic on since the diffusion process is nonsingular. Naturally, the invariant measure, that is, the steady state probability density function of the process , can be determined by the following FPK equation: where and .
The unique solution of (31) satisfies both the periodicity condition and the normality condition , respectively. Thus, the solution of (31) can be expressed as where and is the modified Bessel function of the first kind.
According to Oseledec’s multiplicative ergodic theorem, it can be concluded that, for any initial value, the exponential growth rate (i.e., the Lyapunov exponent) of (29a) and (29b) for the trivial solution can be described as where w.p.1 means with probability one (almost sure).
From (33), we can obtain two different Lyapunov exponents. Thus, the almost sure stability of the trivial response of (33) can be determined by the largest Lyapunov exponent ; that is, is the bifurcation point of the stability of the trivial response. According to [30], the largest Lyapunov exponent is given as
The variation of the largest Lyapunov exponent as surface and the corresponding isohypse curves determined by (34) for the trivial response of the system are shown in Figures 2, 3, and 4 for cases of , 0.1, and 1 and , respectively, where the inner area of the isohypse curves indicates the almost sure instability, whereas the outer area expresses the almost sure stability. From the mesh surfaces in the figures we can find that near the parameter resonance at excitation frequency , is the largest the Lyapunov exponents increase, reaching their maximum values in the center of the instability region. With the increase of , the stability of the trivial response of the system will change. Figures 3(a) and 3(b) correspond to Figures 2(a) and 2(b), respectively. All parameters except the bandwidth of the narrow-band random excitation are kept the same; is increased from 0.01 to 0.1. We can find that Figures 2 and 3 are almost the same; it means that the bandwidth change within 0.01 and 0.1 produces a weak influence on the stability of the system. However, the mesh surface in Figure 4(a) becomes much flatter than that in Figure 3(a) when is 1.0, which implies that the almost sure unstable areas change a lot with the further increase of . In fact, we can find that the bottom of isohypse curve in Figure 4(b) rises compared with that in Figure 3(b), whereas its top is widened, which implies that the increase of may facilitate the almost sure stability of the trivial response and stabilize the system for a lower excitation amplitude but intensify the instability of the trivial response for a higher amplitude . Such results can also be found in Figure 5 for three isohypse curves of resulting from different ; that is, , , and , respectively.
(a)
(b)
(a)
(b)
(a)
(b)
2.2.2. Stochastic Jump and Bifurcation
As the reduced (28b) has the coupled term , the following FPK equation directly derived from (28a) and (28b) will not contain the terms of stationary trivial solution. To circumvent this difficulty and investigate the stochastic jumps and bifurcations between the stationary trivial and nontrivial solutions, normalization method is adopted by introducing the transformation [32] and, introducing it into (28a) and (28b), we havewhere and are the two perpendicular components of ; that is, .
Equations (36a) and (36b) are equivalent to the following set of Itô equation:where and the last terms in and are the Wong-Zakai correction terms.
Equations (37a) and (37b) are a two-dimensional diffusion process and the corresponding FPK equation is described as given below where is the transition probability density of and .
The initial condition of FPK equation (39) is The boundary conditions with respect to and are
When and , (28a) and (28b) is a deterministically excited system and it has one trivial solution and two possible steady state nontrivial solutions when as given below
The theoretical frequency-response curves based on (42) are shown in Figure 6 for the aforementioned beam when and . Similar to [32], there are also three distinct regions in Figure 6; in region I, only the trivial solution is possible and it is stable; in region II, there are one unstable trivial solution and one stable nontrivial solution; and in region III, there are one stable trivial solution, one stable nontrivial solution, and one unstable nontrivial solution.
Figure 7 lists a series of change of the stationary joint probability density of and , which is numerically solved from the FPK equation (39) by using finite difference method [34, 35] to different excitation central frequency when and . Numerical investigation shows that the numerical result becomes convergent when the grid number is greater than within the amplitude region ; thus we choose the grid number as the following numerical grid parameter. From Figure 7 we can find that the joint probability density has two possible types of peaks in region III: an outer flabellate peak and a central volcano peak, which implies an oscillating motion around the trivial branch. The former represents one more probable motion around the stable stationary nontrivial branch and the latter embodies another more probable motion around the stable stationary trivial branch. Concretely, when the excitation central frequency locates in region II, the more probable motion is basically around the stable stationary nontrivial branch (see Figure 7(a)). As the excitation central frequency increases and enters region III, the transition of the response from the stable stationary nontrivial branch to the stable stationary trivial branch occurs (see Figure 7(b)) and gradually becomes stronger (see Figures 7(c), 7(d), and 7(e), resp.). Finally, when the excitation central frequency reaches or exceeds a certain value, the transition process finishes and the more probable motion is around the stable stationary trivial branch (see Figure 7(f)). In fact, the transition of the response in region III between two stable stationary branches occurs again and again, but the trend is that the stochastic jump is towards the stable stationary trivial branch.
| (a) |
| (b) |
| (c) |
| (d) |
| (e) |
| (f) |
Figure 8 reveals the stationary joint probability density to different bandwidths when and . Hear the choice of means that the excitation central frequency lies in region III in Figure 6. We can find that the joint probability density has two possible types of peaks: an outer flabellate peak and a central peak, which is different from the central volcano peak and means an almost dead motion towards the trivial branch. When the bandwidth is a smaller value, for instance, in Figure 8(a), the more probable motion is around the stable stationary nontrivial branch and there is only a sharp and narrow flabellate peak. As the bandwidth increases, the outer flabellate peak gradually becomes a short and wide one (see Figures 8(b) and 8(c)). Finally, when the bandwidth continues increasing, the transition process will finish and the more probable motion will die at the stable stationary trivial branch (see Figure 8(d)).
| (a) |
| (b) |
| (c) |
| (d) |
The theoretical force-response curve based on (42) is shown in Figure 9 to make further research on the effect of excitation force. It can also be found from Figure 9 that there are three distinct regions: region I, that is, the single-valued region, only the trivial solution is possible and it is stable; region II, that is, the triple-valued region, there are one stable trivial solution, one stable nontrivial solution, and one unstable nontrivial solution; and region III, that is, the dual-valued region, there are one unstable trivial solution and one stable nontrivial solution.
Figure 10 shows the stationary joint probability density of the system to different excitation forces when and . We can find that there is only an outer flabellate peak when the excitation force, for instance, , lies in the region III (see Figure 10(a)). When the excitation force has hardly entered the region II, the stationary joint probability density has two peaks: an outer flabellate peak and a weak central volcano peak, which implies that jumps may occur in this region (see Figure 10(b)), but the corresponding parameters, such as bandwidth and force, are not enough to stimulate the transition of the response. However, as the value of the excitation force decreases, the transition of the response from the stable stationary nontrivial branch to the stable stationary trivial branch becomes stronger (see Figures 10(c) and 10(d)). Such phenomena imply that the lower the excitation force is, the more probable the jump from the stable stationary nontrivial branch to the stable stationary trivial one in region II is.
| (a) |
| (b) |
| (c) |
| (d) |
In order to bring forth the influence of the layup or laminate properties, in what follows, the fiber orientation angle of the first lamina with respect to the -axis of the beam is taken into consideration. Table 1 lists a series of data for the dimensionless natural frequency , the dimensionless excitation amplitude , the dimensionless central frequency , and the nonlinear parameter with a smaller change of the fiber orientation angle . In order to make uniform comparison, the excitation parameters are kept the same; in other words, the values of the excitation parameters such as the central frequency , the amplitude , and the bandwidth in (18) are fixed; here we choose , , and .
Figure 11 shows the stationary joint probability density of the system to different fiber orientation angles when , , and . The overall tendency of the response of the system is that the probable motion jumps towards the stable stationary trivial branch as the fiber orientation angle increases. Concretely, when (see Figures 11(a)–11(c)), the distribution of the stationary joint probability density of the system is almost the same; that is, most of the probable motion is around the stable nontrivial solution. When , the probable motion slightly moves towards the stable stationary nontrivial branch (see Figure 11(d)). As the fiber orientation angle increases from to , the probable motion gradually jumps towards the stable stationary trivial branch (see Figures 11(e)–11(j)), whereas, when the fiber orientation angle is in the region of (), the probable motion goes back towards the stable stationary nontrivial branch (see Figures 11(k)-11(l)). When , the probable motion around the stable stationary trivial branch almost disappears (see Figure 11(m)). Then, as the fiber orientation angle increases, the probable motion strongly jumps towards the stable stationary trivial branch (see Figures 11(n)–11(p)). Finally, the most probable motion is around the stable stationary trivial solution when (see Figure 11(q)).
| (a) |
| (b) |
| (c) |
| (d) |
| (e) |
| (f) |
| (g) |
| (h) |
| (i) |
| (j) |
| (k) |
| (l) |
| (m) |
| (n) |
| (o) |
| (p) |
| (q) |
In order to confirm the validity of the approximate analytical solution, comparison has been made between the approximate solution and the direct numerical simulation. Also, in order to make further and uniform comparison, the stationary probability density of amplitude is taken into consideration. Figures 12(a) and 12(b) correspond to Figures 11(f) and 11(p), respectively; that is, all the parameters are kept the same in both figures. The solid lines in both Figures 12(a) and 12(b) are obtained from FPK equation (39) by projecting the amplitude components and in both Figures 11(f) and 11(p) to the amplitude , whereas the solid squares are numerically integrated from (28a) and (28b) by the pseudorandom signal simulation [36, 37]. It can be seen from Figure 12 that the probability densities obtained from the finite difference method are in good agreement with those from digital simulation of the nonlinear modulation equations.
| (a) |
| (b) |
3. Conclusions
Results show that the mesh surface of the largest Lyapunov exponent becomes flatter as increases, which indicates that the increase of may facilitate the almost sure stability of the trivial response and stabilize the system for a lower excitation force. The basic phenomena indicate that the most probable motion is around the stationary nontrivial solution when the bandwidth is smaller, whereas the most probable motion gradually approaches the stationary trivial solution when the bandwidth becomes higher. If the excitation central frequency is a variable, the basic phenomena imply that the higher the excitation central frequency is, the more probable the jump from the stationary nontrivial solution to the stationary trivial one is once the frequency exceeds a certain value. Concretely, as the value of the excitation force decreases, the outer flabellate peak decreases, while the central volcano peak increases.
In order to bring forth the influence of the layup or laminate properties, the fiber orientation angle of the first lamina with respect to the -axis of the beam is taken into consideration. The overall tendency of the response of the system is that the probable motion jumps from the stable stationary nontrivial branch to the stable stationary trivial one as increases from to .
Finally, comparison has been made between the approximate solution and the direct numerical simulation, which shows that the probability densities obtained from the finite difference method are in good agreement with those from the digital simulation of the nonlinear modulation equations.
Appendix
All coefficients obtained in (2a), (2b), and (2c) are presented as follows: where , and denotes the transformed reduced stiffness constant, for the th orthotropic lamina, it can be written aswhere , is the fiber orientation angle of the th lamina with respect to the -axis of the beam, and is the transformed reduced stiffness in the principal material direction, which is expressed as follows: where , , and are the engineering parameters of the th lamina.
Conflict of Interests
The authors declare that there is no conflict of interests regarding the publication of this paper.
Acknowledgments
The authors gratefully acknowledge the support of the National Natural Science Foundation of China (NNSFC) through Grant no. 11072164 and a Project Funded by the Priority Academic Program Development of Jiangsu Higher Education Institutions (PAPD).