Stochastic Principal Parametric Resonances of Composite Laminated Beams

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 x -axis of the beam increases from zero to a smaller angle.

Xiangjun Lan 1 and Zhihua Feng 1 and Fan Lv 2

Academic Editor:Ahmet S. Yigit

1, School of Mechanical and Electric Engineering, Soochow University, Suzhou 215021, China
2, School of Shagang Iron and Steel, Soochow University, Suzhou 215021, China

Received 17 December 2013; Revised 6 March 2014; Accepted 6 March 2014; 14 April 2014

This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

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 [figure omitted; refer to PDF] and thickness [figure omitted; refer to PDF] subject to axial load [figure omitted; refer to PDF] . The displacement proposed by Levinson [4] and used by Reddy [6] is given as below [figure omitted; refer to PDF] where [figure omitted; refer to PDF] , [figure omitted; refer to PDF] and [figure omitted; refer to PDF] are total displacements, [figure omitted; refer to PDF] and [figure omitted; refer to PDF] are the displacements of a point in the middle surface of the beam in the [figure omitted; refer to PDF] and [figure omitted; refer to PDF] directions, and [figure omitted; refer to PDF] is rotation of the normal about the [figure omitted; refer to PDF] -axis.

Figure 1: Configuration of a laminated composite beam.

[figure omitted; refer to PDF]

According to Reddy and von Karman's theory, the governing equations of the motion by the Hamilton principle can be obtained as [figure omitted; refer to PDF] where [figure omitted; refer to PDF] is the coefficient of damping and other coefficients can be found in the appendix.

Also, the boundary conditions can be written as

: [figure omitted; refer to PDF] : [figure omitted; refer to PDF]

: [figure omitted; refer to PDF] : [figure omitted; refer to PDF]

: [figure omitted; refer to PDF] and [figure omitted; refer to PDF] : [figure omitted; refer to PDF]

: [figure omitted; refer to PDF] and [figure omitted; refer to PDF] : [figure omitted; refer to PDF]

: [figure omitted; refer to PDF] and [figure omitted; refer to PDF] : [figure omitted; refer to PDF]

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 [figure omitted; refer to PDF]

Substituting (7) into (2a) and (2c) and neglecting the longitudinal inertia and the moment of inertia about the neutral axis, we have [figure omitted; refer to PDF]

Integrating (8) and using (3a) and (3b), we arrive at [figure omitted; refer to PDF]

According to the boundary conditions, the solution of (9) can be approximately expressed as [33] [figure omitted; refer to PDF]

In order to obtain the dimensionless equations, the variables and parameters are introduced as given below [figure omitted; refer to PDF]

Naturally, for simplicity, the star is dropped in the following analysis.

Closed form solutions for [figure omitted; refer to PDF] can be obtained for the beam by assuming [figure omitted; refer to PDF]

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 [figure omitted; refer to PDF] 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 [figure omitted; refer to PDF] where [figure omitted; refer to PDF]

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, [figure omitted; refer to PDF] where [figure omitted; refer to PDF] and [figure omitted; refer to PDF] are the amplitude and center frequency of the excitation, which are constants, [figure omitted; refer to PDF] is an intensity, which represents the bandwidth, of the narrow-band random excitation, [figure omitted; refer to PDF] is a standard Wiener process, and [figure omitted; refer to PDF] is a uniformly distributed random number in [figure omitted; refer to PDF] . The inclusion of the phase angle [figure omitted; refer to PDF] makes [figure omitted; refer to PDF] a stationary process. So [figure omitted; refer to PDF] is a non-Gaussian distributed stationary stochastic process and has density function [figure omitted; refer to PDF] , zero mean value, and spectral density function as follows: [figure omitted; refer to PDF]

The bandwidth of process [figure omitted; refer to PDF] depends mainly on parameter [figure omitted; refer to PDF] . It is a narrow-band process when [figure omitted; refer to PDF] is small and a wide-band process when [figure omitted; refer to PDF] is large. Thus, (14) can be rewritten as [figure omitted; refer to PDF] where [figure omitted; refer to PDF] , [figure omitted; refer to PDF] , [figure omitted; refer to PDF] , and [figure omitted; refer to PDF] is a small parameter.

Also, a new time scale is defined as [figure omitted; refer to PDF] . Let [figure omitted; refer to PDF] , [figure omitted; refer to PDF] , [figure omitted; refer to PDF] , and [figure omitted; refer to PDF] and rewrite (18) as [figure omitted; refer to PDF]

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 [figure omitted; refer to PDF] m, [figure omitted; refer to PDF] m, [figure omitted; refer to PDF] m, [figure omitted; refer to PDF] kg/m³ , [figure omitted; refer to PDF] GPa, [figure omitted; refer to PDF] GPa, [figure omitted; refer to PDF] GPa, [figure omitted; refer to PDF] GPa, [figure omitted; refer to PDF] , and [figure omitted; refer to PDF] . Thus, parameter [figure omitted; refer to PDF] in (19) is calculated to be 2.4566 when the fiber orientation angle [figure omitted; refer to PDF] of the first lamina with respect to the [figure omitted; refer to PDF] -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 [figure omitted; refer to PDF]

Moreover, for unit Wiener progress [figure omitted; refer to PDF] , because of [figure omitted; refer to PDF] and [figure omitted; refer to PDF] , we have [figure omitted; refer to PDF]

Here the discussion and investigation are restricted to the case of principal parametric resonance of the first mode. Thus, the frequency detuning parameter [figure omitted; refer to PDF] is introduced to describe the closeness to the principal parametric resonance as given below [figure omitted; refer to PDF] Substituting (20) and (22) into (19) yields

: order [figure omitted; refer to PDF] : [figure omitted; refer to PDF]

: order [figure omitted; refer to PDF] : [figure omitted; refer to PDF]

The first-order approximate solution of (23) is [figure omitted; refer to PDF] where cc stands for the complex conjugate of the preceding terms.

Substituting (22) and (25) into (24) and eliminating the secular terms, we have [figure omitted; refer to PDF] where [figure omitted; refer to PDF] .

Substituting the polar form [figure omitted; refer to PDF] into (26) and separating the real and imaginary parts, we get [figure omitted; refer to PDF]

2.2.1. Stability Analysis of the Trivial Response

It can be found that (28a) and (28b) have a trivial steady state solution [figure omitted; refer to PDF] . In order to determine the stability of the trivial steady state response, the trivial response is expanded near [figure omitted; refer to PDF] and the corresponding nonlinear terms are neglected; the linearization form of (28a) and (28b) can finally be expressed as [figure omitted; refer to PDF]

Let [figure omitted; refer to PDF] and we rewrite (29a) and (29b) in the form of Itô equation as given below [figure omitted; refer to PDF]

It can be found that the stochastic process [figure omitted; refer to PDF] generated on [figure omitted; refer to PDF] by (30a) and (30b) is Markov and is also ergodic on [figure omitted; refer to PDF] since the diffusion process is nonsingular. Naturally, the invariant measure, that is, the steady state probability density function [figure omitted; refer to PDF] of the process [figure omitted; refer to PDF] , can be determined by the following FPK equation: [figure omitted; refer to PDF] where [figure omitted; refer to PDF] and [figure omitted; refer to PDF] .

The unique solution of (31) satisfies both the periodicity condition [figure omitted; refer to PDF] and the normality condition [figure omitted; refer to PDF] , respectively. Thus, the solution of (31) can be expressed as [figure omitted; refer to PDF] where [figure omitted; refer to PDF] and [figure omitted; refer to PDF] 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 [figure omitted; refer to PDF] 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 [figure omitted; refer to PDF] ; that is, [figure omitted; refer to PDF] is the bifurcation point of the stability of the trivial response. According to [30], the largest Lyapunov exponent [figure omitted; refer to PDF] is given as [figure omitted; refer to PDF]

The variation of the largest Lyapunov exponent [figure omitted; refer to PDF] as [figure omitted; refer to PDF] 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 [figure omitted; refer to PDF] , 0.1, and 1 and [figure omitted; refer to PDF] , 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 [figure omitted; refer to PDF] , is the largest the Lyapunov exponents increase, reaching their maximum values in the center of the instability region. With the increase of [figure omitted; refer to PDF] , 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; [figure omitted; refer to PDF] 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 [figure omitted; refer to PDF] is 1.0, which implies that the almost sure unstable areas change a lot with the further increase of [figure omitted; refer to PDF] . In fact, we can find that the bottom of isohypse curve [figure omitted; refer to PDF] in Figure 4(b) rises compared with that in Figure 3(b), whereas its top is widened, which implies that the increase of [figure omitted; refer to PDF] may facilitate the almost sure stability of the trivial response and stabilize the system for a lower excitation amplitude [figure omitted; refer to PDF] but intensify the instability of the trivial response for a higher amplitude [figure omitted; refer to PDF] . Such results can also be found in Figure 5 for three isohypse curves of [figure omitted; refer to PDF] resulting from different [figure omitted; refer to PDF] ; that is, [figure omitted; refer to PDF] , [figure omitted; refer to PDF] , and [figure omitted; refer to PDF] , respectively.

The largest Lyapunov exponent [figure omitted; refer to PDF] of the trivial response of the system: [figure omitted; refer to PDF] , [figure omitted; refer to PDF] . (a) Mesh surface of [figure omitted; refer to PDF] ; (b) isohypse curves of [figure omitted; refer to PDF] .

(a) [figure omitted; refer to PDF]

(b) [figure omitted; refer to PDF]

The largest Lyapunov exponent [figure omitted; refer to PDF] of the trivial response of the system: [figure omitted; refer to PDF] , [figure omitted; refer to PDF] . (a) Mesh surface of [figure omitted; refer to PDF] ; (b) isohypse curves of [figure omitted; refer to PDF] .

(a) [figure omitted; refer to PDF]

(b) [figure omitted; refer to PDF]

The largest Lyapunov exponent [figure omitted; refer to PDF] of the trivial response of the system: [figure omitted; refer to PDF] , [figure omitted; refer to PDF] . (a) Mesh surface of [figure omitted; refer to PDF] ; (b) isohypse curves of [figure omitted; refer to PDF] .

(a) [figure omitted; refer to PDF]

(b) [figure omitted; refer to PDF]

Figure 5: Isohypse curves of [figure omitted; refer to PDF] .

[figure omitted; refer to PDF]

2.2.2. Stochastic Jump and Bifurcation

As the reduced (28b) has the coupled term [figure omitted; refer to PDF] , 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] [figure omitted; refer to PDF] and, introducing it into (28a) and (28b), we have [figure omitted; refer to PDF] where [figure omitted; refer to PDF] and [figure omitted; refer to PDF] are the two perpendicular components of [figure omitted; refer to PDF] ; that is, [figure omitted; refer to PDF] .

Equations (36a) and (36b) are equivalent to the following set of Itô equation: [figure omitted; refer to PDF] where [figure omitted; refer to PDF] and the last terms in [figure omitted; refer to PDF] and [figure omitted; refer to PDF] 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 [figure omitted; refer to PDF] where [figure omitted; refer to PDF] is the transition probability density of [figure omitted; refer to PDF] and [figure omitted; refer to PDF] .

The initial condition of FPK equation (39) is [figure omitted; refer to PDF] The boundary conditions with respect to [figure omitted; refer to PDF] and [figure omitted; refer to PDF] are [figure omitted; refer to PDF]

When [figure omitted; refer to PDF] and [figure omitted; refer to PDF] , (28a) and (28b) is a deterministically excited system and it has one trivial solution and two possible steady state nontrivial solutions when [figure omitted; refer to PDF] as given below [figure omitted; refer to PDF]

The theoretical frequency-response curves based on (42) are shown in Figure 6 for the aforementioned beam when [figure omitted; refer to PDF] and [figure omitted; refer to PDF] . 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 6: Frequency-response curves when [figure omitted; refer to PDF] and [figure omitted; refer to PDF] .

[figure omitted; refer to PDF]

Figure 7 lists a series of change of the stationary joint probability density of [figure omitted; refer to PDF] and [figure omitted; refer to PDF] , which is numerically solved from the FPK equation (39) by using finite difference method [34, 35] to different excitation central frequency [figure omitted; refer to PDF] when [figure omitted; refer to PDF] and [figure omitted; refer to PDF] . Numerical investigation shows that the numerical result becomes convergent when the grid number is greater than [figure omitted; refer to PDF] within the amplitude region [figure omitted; refer to PDF] ; thus we choose the grid number [figure omitted; refer to PDF] 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 [figure omitted; refer to PDF] 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 [figure omitted; refer to PDF] 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 [figure omitted; refer to PDF] 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.

Stationary joint probability density of the system to different excitation central frequencies [figure omitted; refer to PDF] : [figure omitted; refer to PDF] , [figure omitted; refer to PDF] . (a) [figure omitted; refer to PDF] ; (b) [figure omitted; refer to PDF] ; (c) [figure omitted; refer to PDF] ; (d) [figure omitted; refer to PDF] ; (e) [figure omitted; refer to PDF] ; (f) [figure omitted; refer to PDF] .

(a) [figure omitted; refer to PDF]

[figure omitted; refer to PDF]

(b) [figure omitted; refer to PDF]

[figure omitted; refer to PDF]

(c) [figure omitted; refer to PDF]

[figure omitted; refer to PDF]

(d) [figure omitted; refer to PDF]

[figure omitted; refer to PDF]

(e) [figure omitted; refer to PDF]

[figure omitted; refer to PDF]

(f) [figure omitted; refer to PDF]

[figure omitted; refer to PDF]

Figure 8 reveals the stationary joint probability density to different bandwidths when [figure omitted; refer to PDF] and [figure omitted; refer to PDF] . Hear the choice of [figure omitted; refer to PDF] 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 [figure omitted; refer to PDF] is a smaller value, for instance, [figure omitted; refer to PDF] 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 [figure omitted; refer to PDF] increases, the outer flabellate peak gradually becomes a short and wide one (see Figures 8(b) and 8(c)). Finally, when the bandwidth [figure omitted; refer to PDF] continues increasing, the transition process will finish and the more probable motion will die at the stable stationary trivial branch (see Figure 8(d)).

Stationary joint probability density of the system to different bandwidths [figure omitted; refer to PDF] . [figure omitted; refer to PDF] , [figure omitted; refer to PDF] .

(a) [figure omitted; refer to PDF]

[figure omitted; refer to PDF]

(b) [figure omitted; refer to PDF]

[figure omitted; refer to PDF]

(c) [figure omitted; refer to PDF]

[figure omitted; refer to PDF]

(d) [figure omitted; refer to PDF]

[figure omitted; refer to PDF]

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 9: Force-response curves when [figure omitted; refer to PDF] .

[figure omitted; refer to PDF]

Figure 10 shows the stationary joint probability density of the system to different excitation forces when [figure omitted; refer to PDF] and [figure omitted; refer to PDF] . We can find that there is only an outer flabellate peak when the excitation force, for instance, [figure omitted; refer to PDF] , 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.

Stationary joint probability density of the system to different excitation forces. [figure omitted; refer to PDF] , [figure omitted; refer to PDF] .

(a) [figure omitted; refer to PDF]

[figure omitted; refer to PDF]

(b) [figure omitted; refer to PDF]

[figure omitted; refer to PDF]

(c) [figure omitted; refer to PDF]

[figure omitted; refer to PDF]

(d) [figure omitted; refer to PDF]

[figure omitted; refer to PDF]

In order to bring forth the influence of the layup or laminate properties, in what follows, the fiber orientation angle [figure omitted; refer to PDF] of the first lamina with respect to the [figure omitted; refer to PDF] -axis of the beam is taken into consideration. Table 1 lists a series of data for the dimensionless natural frequency [figure omitted; refer to PDF] , the dimensionless excitation amplitude [figure omitted; refer to PDF] , the dimensionless central frequency [figure omitted; refer to PDF] , and the nonlinear parameter [figure omitted; refer to PDF] with a smaller change of the fiber orientation angle [figure omitted; refer to PDF] . 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 [figure omitted; refer to PDF] , the amplitude [figure omitted; refer to PDF] , and the bandwidth [figure omitted; refer to PDF] in (18) are fixed; here we choose [figure omitted; refer to PDF] , [figure omitted; refer to PDF] , and [figure omitted; refer to PDF] .

Table 1: System and excitation parameters corresponding to [figure omitted; refer to PDF] .

[figure omitted; refer to PDF]	[figure omitted; refer to PDF]	[figure omitted; refer to PDF]	[figure omitted; refer to PDF]	[figure omitted; refer to PDF]	B
0°	1.8067	8.0191	2.4566	1.0000	3.0000
1°	1.8062	8.0144	2.4566	1.0006	3.0008
2°	1.8046	8.0004	2.4566	1.0023	3.0035
3°	1.8020	7.9771	2.4566	1.0052	3.0078
4°	1.7983	7.9445	2.4566	1.0094	3.0140
5°	1.7936	7.9028	2.4565	1.0147	3.0220
6°	1.7879	7.8522	2.4565	1.0211	3.0316
7°	1.7811	7.7927	2.4564	1.0290	3.0431
8°	1.7733	7.7246	2.4563	1.0380	3.0565
9°	1.7646	7.6480	2.4562	1.0483	3.0715
10°	1.7548	7.5634	2.4561	1.0600	3.0887
11°	1.7441	7.4708	2.4559	1.0731	3.1076
12°	1.7324	7.3707	2.4558	1.0876	3.1286
13°	1.7199	7.2633	2.4555	1.1035	3.1514
14°	1.7064	7.1490	2.4553	1.1210	3.1763
15°	1.6920	7.0281	2.4549	1.1402	3.2034
16°	1.6768	6.9011	2.4546	1.1609	3.2324

Figure 11 shows the stationary joint probability density of the system to different fiber orientation angles when [figure omitted; refer to PDF] , [figure omitted; refer to PDF] , and [figure omitted; refer to PDF] . 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 [figure omitted; refer to PDF] (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 [figure omitted; refer to PDF] , the probable motion slightly moves towards the stable stationary nontrivial branch (see Figure 11(d)). As the fiber orientation angle increases from [figure omitted; refer to PDF] to [figure omitted; refer to PDF] , 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 ( [figure omitted; refer to PDF] ), the probable motion goes back towards the stable stationary nontrivial branch (see Figures 11(k)-11(l)). When [figure omitted; refer to PDF] , 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 [figure omitted; refer to PDF] (see Figure 11(q)).

Stationary joint probability density of the system to different fiber orientation angles [figure omitted; refer to PDF] of the first lamina with respect to the [figure omitted; refer to PDF] -axis of the beam. [figure omitted; refer to PDF] , [figure omitted; refer to PDF] , and [figure omitted; refer to PDF] .

(a) [figure omitted; refer to PDF]

[figure omitted; refer to PDF]

(b) [figure omitted; refer to PDF]

[figure omitted; refer to PDF]

(c) [figure omitted; refer to PDF]

[figure omitted; refer to PDF]

(d) [figure omitted; refer to PDF]

[figure omitted; refer to PDF]

(e) [figure omitted; refer to PDF]

[figure omitted; refer to PDF]

(f) [figure omitted; refer to PDF]

[figure omitted; refer to PDF]

(g) [figure omitted; refer to PDF]

[figure omitted; refer to PDF]

(h) [figure omitted; refer to PDF]

[figure omitted; refer to PDF]

(i) [figure omitted; refer to PDF]

[figure omitted; refer to PDF]

(j) [figure omitted; refer to PDF]

[figure omitted; refer to PDF]

(k) [figure omitted; refer to PDF]

[figure omitted; refer to PDF]

(l) [figure omitted; refer to PDF]

[figure omitted; refer to PDF]

(m) [figure omitted; refer to PDF]

[figure omitted; refer to PDF]

(n) [figure omitted; refer to PDF]

[figure omitted; refer to PDF]

(o) [figure omitted; refer to PDF]

[figure omitted; refer to PDF]

(p) [figure omitted; refer to PDF]

[figure omitted; refer to PDF]

(q) [figure omitted; refer to PDF]

[figure omitted; refer to PDF]

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 [figure omitted; refer to PDF] and [figure omitted; refer to PDF] in both Figures 11(f) and 11(p) to the amplitude [figure omitted; refer to PDF] , 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.

Stationary probability density of the amplitude to different fiber orientation angles [figure omitted; refer to PDF] of the first lamina with respect to the [figure omitted; refer to PDF] -axis of the beam, [figure omitted; refer to PDF] , [figure omitted; refer to PDF] , and [figure omitted; refer to PDF] .

(a) [figure omitted; refer to PDF]

[figure omitted; refer to PDF]

(b) [figure omitted; refer to PDF]

[figure omitted; refer to PDF]

3. Conclusions

Results show that the mesh surface of the largest Lyapunov exponent becomes flatter as [figure omitted; refer to PDF] increases, which indicates that the increase of [figure omitted; refer to PDF] 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 [figure omitted; refer to PDF] of the first lamina with respect to the [figure omitted; refer to PDF] -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 [figure omitted; refer to PDF] increases from [figure omitted; refer to PDF] to [figure omitted; refer to PDF] .

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.

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).

Conflict of Interests

The authors declare that there is no conflict of interests regarding the publication of this paper.

References

[1] O. C. Zinenkiewicz, R. L. Taylor The Finite Element Method , vol. 1, McGraw-Hill, Singapore, 1989.

[2] E. Reissner, "On transverse bending of plates, including the effect of transverse shear deformation," International Journal of Solids and Structures , vol. 11, no. 5, pp. 569-573, 1975.

[3] R. D. Mindlin, "Influence of rotatory inertia and shear on flexural motions of isotropic, elastic plates," Journal of Applied Mechanics , vol. 18, pp. 31-38, 1951.

[4] M. Levinson, "A new rectangular beam theory," Journal of Sound and Vibration , vol. 74, no. 1, pp. 81-87, 1981.

[5] A. V. K. Murty, "Toward a consistent beam theory," AIAA Journal , vol. 22, no. 6, pp. 811-816, 1984.

[6] J. N. Reddy, "A simple higher-order theory for laminated composite plates," Journal of Applied Mechanics, Transactions ASME , vol. 51, no. 4, pp. 745-752, 1984.

[7] T. Kant, S. R. Marur, G. S. Rao, "Analytical solution to the dynamic analysis of laminated beams using higher order refined theory," Composite Structures , vol. 40, no. 1, pp. 1-9, 1997.

[8] H. Matsunaga, "Vibration and buckling of multilayered composite beams according to higher order deformation theories," Journal of Sound and Vibration , vol. 246, no. 1, pp. 47-62, 2001.

[9] M. Touratier, "An efficient standard plate theory," International Journal of Engineering Science , vol. 29, no. 8, pp. 901-916, 1991.

[10] M. Touratier, "A refined theory of laminated shallow shells," International Journal of Solids and Structures , vol. 29, no. 11, pp. 1401-1415, 1992.

[11] P. Vidal, O. Polit, "A family of sinus finite elements for the analysis of rectangular laminated beams," Composite Structures , vol. 84, no. 1, pp. 56-72, 2008.

[12] M. Karama, K. S. Afaq, S. Mistou, "Mechanical behaviour of laminated composite beam by the new multi-layered laminated composite structures model with transverse shear stress continuity," International Journal of Solids and Structures , vol. 40, no. 6, pp. 1525-1546, 2003.

[13] M. Aydogdu, "A new shear deformation theory for laminated composite plates," Composite Structures , vol. 89, no. 1, pp. 94-101, 2009.

[14] S. J. Song, A. M. Waas, "Effects of shear deformation on buckling and free vibration of laminated composite beams," Composite Structures , vol. 37, no. 1, pp. 33-43, 1997.

[15] P. Subramanian, "Dynamic analysis of laminated composite beams using higher order theories and finite elements," Composite Structures , vol. 73, no. 3, pp. 342-353, 2006.

[16] M. R. M. Crespo Da Silva, C. C. Glynn, "Non-linear flexural-flexural-torsional dynamics of inextensional beams-I: equations of motion," Journal of Structural Mechanics , vol. 6, no. 4, pp. 437-448, 1978.

[17] M. R. M. Crespo Da Silva, C. C. Glynn, "Non-linear flexural-flexural-torsional dynamics of inextensional beams-II: forced motions," Journal of Structural Mechanics , vol. 6, no. 4, pp. 449-461, 1978.

[18] A. H. Nayfeh, P. F. Pai, "Non-linear non-planar parametric responses of an inextensional beam," International Journal of Non-Linear Mechanics , vol. 24, no. 2, pp. 139-158, 1989.

[19] P. F. Pai, A. H. Nayfeh, "Non-linear non-planar oscillations of a cantilever beam under lateral base excitations," International Journal of Non-Linear Mechanics , vol. 25, no. 5, pp. 455-474, 1990.

[20] L. D. Zavodney, A. H. Nayfeh, "The non-linear response of a slender beam carrying a lumped mass to a principal parametric excitation: theory and experiment," International Journal of Non-Linear Mechanics , vol. 24, no. 2, pp. 105-125, 1989.

[21] R. C. Kar, S. K. Dwivedy, "Non-linear dynamics of a slender beam carrying a lumped mass with principal parametric and internal resonances," International Journal of Non-Linear Mechanics , vol. 34, no. 3, pp. 515-529, 1998.

[22] S. K. Dwivedy, R. C. Kar, "Dynamics of a slender beam with an attached mass under combination parametric and internal resonances part I: steady state response," Journal of Sound and Vibration , vol. 221, no. 5, pp. 823-848, 1999.

[23] S. K. Dwivedy, R. C. Kar, "Dynamics of a slender beam with an attached mass under combination parametric and internal resonances, part II: periodic and chaotic responses," Journal of Sound and Vibration , vol. 222, no. 2, pp. 281-305, 1999.

[24] S. K. Dwivedy, R. C. Kar, "Non-linear dynamics of a slender beam carrying a lumped mass under principal parametric resonance with three-mode interactions," International Journal of Non-Linear Mechanics , vol. 36, no. 6, pp. 927-945, 2001.

[25] S. K. Dwivedy, R. C. Kar, "Simultaneous combination, principal parametric and internal resonances in a slender beam with a lumped mass: three-mode interactions," Journal of Sound and Vibration , vol. 242, no. 1, pp. 27-46, 2001.

[26] S. K. Dwivedy, R. C. Kar, "Simultaneous combination and 1:3:5 internal resonances in a parametrically excited beam-mass system," International Journal of Non-Linear Mechanics , vol. 38, no. 4, pp. 585-596, 2003.

[27] T. J. Anderson, B. Balachandran, A. H. Nayfeh, "Nonlinear resonances in a flexible cantilever beam," Journal of Vibration and Acoustics, Transactions of the ASME , vol. 116, no. 4, pp. 480-484, 1994.

[28] T. J. Anderson, A. H. Nayfeh, B. Balachandran, "Coupling between high-frequency modes and a low-frequency mode: theory and experiment," Nonlinear Dynamics , vol. 11, no. 1, pp. 17-36, 1996.

[29] T. J. Anderson, A. H. Nayfeh, B. Balachandran, "Experimental verification of the importance of the nonlinear curvature in the response of a cantilever beam," Journal of Vibration and Acoustics, Transactions of the ASME , vol. 118, no. 1, pp. 21-27, 1996.

[30] Z. H. Feng, X. J. Lan, X. D. Zhu, "Principal parametric resonances of a slender cantilever beam subject to axial narrow-band random excitation of its base," International Journal of Non-Linear Mechanics , vol. 42, no. 10, pp. 1170-1185, 2007.

[31] Z. H. Feng, X. J. Lan, X. D. Zhu, "Explanation on the importance of narrow-band random excitation characters in the response of a cantilever beam," Journal of Sound and Vibration , vol. 325, no. 4-5, pp. 923-937, 2009.

[32] Z. H. Feng, X. D. Zhu, X. J. Lan, "Stochastic jump and bifurcation of a slender cantilever beam carrying a lumped mass under narrow-band principal parametric excitation," International Journal of Non-Linear Mechanics , vol. 46, no. 10, pp. 1330-1340, 2011.

[33] Z. Yao, W. Zhang, L. Chen, "Bifurcation, chaotic dynamics and control of piezoelectric laminated composite beam," Chinese Journal of Theoretical and Applied Mechanics , vol. 41, no. 1, pp. 129-140, 2009.

[34] Z. L. Huang, W. Q. Zhu, "Stochastic averaging of quasi-integrable Hamiltonian systems under combined harmonic and white noise excitations," International Journal of Non-Linear Mechanics , vol. 39, no. 9, pp. 1421-1434, 2004.

[35] Z. L. Huang, W. Q. Zhu, "Stochastic averaging of quasi-integrable Hamiltonian systems under bounded noise excitations," Probabilistic Engineering Mechanics , vol. 19, no. 3, pp. 219-228, 2004.

[36] M. Shinozuka, "Simulation of multivariate and multidimensional random processes," Journal of the Acoustical Society of America , vol. 49, no. 1, pp. 357-367, 1971.

[37] M. Shinozuka, C.-M. Jan, "Digital simulation of random processes and its applications," Journal of Sound and Vibration , vol. 25, no. 1, pp. 111-128, 1972.

Appendix

Appendix

All coefficients obtained in (2a), (2b), and (2c) are presented as follows: [figure omitted; refer to PDF] where [figure omitted; refer to PDF] , and [figure omitted; refer to PDF] denotes the transformed reduced stiffness constant, for the [figure omitted; refer to PDF] th orthotropic lamina, it can be written as [figure omitted; refer to PDF] where [figure omitted; refer to PDF] , [figure omitted; refer to PDF] is the fiber orientation angle of the [figure omitted; refer to PDF] th lamina with respect to the [figure omitted; refer to PDF] -axis of the beam, and [figure omitted; refer to PDF] is the transformed reduced stiffness in the principal material direction, which is expressed as follows: [figure omitted; refer to PDF] where [figure omitted; refer to PDF] , [figure omitted; refer to PDF] , and [figure omitted; refer to PDF] are the engineering parameters of the [figure omitted; refer to PDF] th lamina.