# Stochastic Principal Parametric Resonances of Composite Laminated Beams

Authors: Xiangjun Lan; Zhihua Feng; Fan Lv
Journal: Shock and Vibration
 (2014)

---

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


---

## Body

## 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 Figure1 with length L and thickness h subject to axial load f. The displacement proposed by Levinson [4] and used by Reddy [6] is given as below(1a)
u
(
x
,
z
,
t
)
=
u
0
(
x
,
t
)
+
z
φ
x
(
x
,
t
)
-
z
3
c
1
(
φ
x
+
∂
w
0
∂
x
)
,
(1b)
w
(
x
,
z
,
t
)
=
w
0
(
x
,
t
)
,
(1c)
φ
x
=
(
∂
u
∂
z
)
z
=
0
,where c
1
=
4
/
3
h
2, u and w are total displacements, u
0 and w
0 are the displacements of a point in the middle surface of the beam in the x and z directions, and φ
x is rotation of the normal about the y-axis.Figure 1
Configuration of a laminated composite beam.According to Reddy and von Karman’s theory, the governing equations of the motion by the Hamilton principle can be obtained as(2a)
A
11
∂
2
u
0
∂
x
2
+
A
11
∂
w
0
∂
x
∂
2
w
0
∂
x
2
+
(
B
11
-
c
1
E
11
)
∂
2
φ
x
∂
x
2
-
c
1
E
11
∂
3
w
0
∂
x
3
=
I
0
∂
2
u
0
∂
t
2
+
J
1
∂
2
φ
x
∂
t
2
-
c
1
I
3
∂
3
w
0
∂
x
∂
t
2
,
(2b)
A
11
∂
u
0
∂
x
∂
2
w
0
∂
x
2
+
A
11
∂
2
u
0
∂
x
2
∂
w
0
∂
x
+
c
1
E
11
∂
3
u
0
∂
x
3
+
(
A
55
-
6
c
1
D
55
+
9
c
1
2
F
55
)
∂
φ
x
∂
x
+
3
2
A
11
(
∂
w
0
∂
x
)
2
∂
2
w
0
∂
x
2
+
(
B
11
-
c
1
E
11
)
∂
2
φ
x
∂
x
2
∂
w
0
∂
x
+
c
1
(
F
11
-
c
1
H
11
)
∂
3
φ
x
∂
x
3
+
(
A
55
-
6
c
1
D
55
+
9
c
1
2
F
55
)
∂
2
w
0
∂
x
2
+
(
B
11
-
c
1
E
11
)
×
∂
φ
x
∂
x
∂
2
w
0
∂
x
2
-
c
1
2
H
11
∂
4
w
0
∂
x
4
-
c
3
w
˙
0
=
I
0
∂
2
w
0
∂
t
2
+
c
1
I
3
∂
3
u
0
∂
x
∂
t
2
+
c
1
J
4
∂
3
φ
x
∂
x
∂
t
2
-
c
1
2
I
6
∂
4
w
0
∂
x
2
∂
t
2
,
(2c)
(
B
11
-
c
1
E
11
)
∂
2
u
0
∂
x
2
+
(
B
11
-
c
1
E
11
)
∂
w
0
∂
x
∂
2
w
0
∂
x
2
-
(
c
1
F
11
-
c
1
2
H
11
)
∂
3
w
0
∂
x
3
+
[
(
D
11
-
c
1
F
11
)
-
c
1
(
F
11
-
c
1
H
11
)
]
∂
2
φ
x
∂
x
2
-
(
A
55
-
6
c
1
D
55
+
9
c
1
2
F
55
)
φ
x
-
(
A
55
-
6
c
1
D
55
+
9
c
1
2
F
55
)
∂
w
0
∂
x
=
J
1
∂
2
u
0
∂
t
2
+
K
2
∂
2
φ
x
∂
t
2
-
c
1
J
4
∂
3
w
0
∂
x
∂
t
2
,where c
3 is the coefficient of damping and other coefficients can be found in the appendix.Also, the boundary conditions can be written asx
=
0:(3a)
δ
u
0
=
0
or
b
(
A
11
∂
u
0
∂
x
+
1
2
A
11
(
∂
w
0
∂
x
)
2
h
h
h
(
∂
w
0
∂
x
)
2
+
(
B
11
-
c
1
E
11
)
∂
φ
x
∂
x
-
c
1
E
11
∂
2
w
0
∂
x
2
)
+
f
h
=
0
,
x
=
L:
(3b)
δ
u
0
=
0
or
A
11
∂
u
0
∂
x
+
1
2
A
11
(
∂
w
0
∂
x
)
2
+
(
B
11
-
c
1
E
11
)
∂
φ
x
∂
x
-
c
1
E
11
∂
2
w
0
∂
x
2
=
0
,
x
=
0 and x
=
L:(4)
δ
φ
x
=
0
or
(
B
11
z
-
c
1
E
11
z
3
)
[
∂
u
0
∂
x
+
1
2
(
∂
w
0
∂
x
)
2
]
+
(
D
11
-
2
c
1
F
11
+
c
1
2
H
11
)
∂
φ
x
∂
x
-
(
c
1
F
11
-
c
1
2
H
11
)
∂
2
w
0
∂
x
2
=
0
,
x
=
0 and x
=
L:
(5)
∂
δ
w
0
∂
x
=
0
or
E
11
∂
u
0
∂
x
+
1
2
E
11
(
∂
w
0
∂
x
)
2
+
c
1
(
F
11
-
c
1
H
11
)
∂
φ
x
∂
x
-
c
1
2
H
11
∂
2
w
0
∂
x
2
=
0
,
x
=
0 and x
=
L:
(6)
δ
w
0
=
0
,
or
[
A
11
∂
u
0
∂
x
+
1
2
A
11
(
∂
w
0
∂
x
)
2
h
.
+
(
B
11
-
c
1
E
11
)
∂
φ
x
∂
x
-
c
1
E
11
∂
2
w
0
∂
x
2
]
∂
w
0
∂
x
+
c
1
[
E
11
∂
3
u
0
∂
x
3
+
E
11
(
∂
2
w
0
∂
x
2
)
2
+
E
11
∂
w
0
∂
x
∂
3
w
0
∂
x
3
h
.
h
.
h
.
h
.
+
(
F
11
-
c
1
H
11
)
∂
3
φ
x
∂
x
3
-
c
1
H
11
∂
4
w
0
∂
x
4
]
+
(
A
55
-
6
c
1
D
55
+
9
c
1
2
F
55
)
(
φ
x
+
∂
w
0
∂
x
)
+
c
1
I
3
∂
2
u
0
∂
t
2
-
c
1
(
I
4
-
c
1
I
6
)
∂
2
φ
x
∂
t
2
+
c
1
2
I
6
∂
3
w
0
∂
x
∂
t
2
=
0
.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(7)
B
11
=
E
11
=
0
,
I
1
=
I
3
=
0
.Substituting (7) into (2a) and (2c) and neglecting the longitudinal inertia and the moment of inertia about the neutral axis, we have
(8)
A
11
∂
2
u
0
∂
x
2
+
A
11
∂
w
0
∂
x
∂
2
w
0
∂
x
2
=
0
,
(9)
-
(
A
55
-
6
c
1
D
55
+
9
c
1
2
F
55
)
∂
w
0
∂
x
-
(
A
55
-
6
c
1
D
55
+
9
c
1
2
F
55
)
φ
x
+
[
(
D
11
-
c
1
F
11
)
-
c
1
(
F
11
-
c
1
H
11
)
]
∂
2
φ
x
∂
x
2
-
(
c
1
F
11
-
c
1
2
H
11
)
∂
3
w
0
∂
x
3
=
0
.Integrating (8) and using (3a) and (3b), we arrive at
(10)
∂
u
0
∂
x
=
-
1
2
(
∂
w
0
∂
x
)
2
+
h
f
A
11
.According to the boundary conditions, the solution of (9) can be approximately expressed as [33]
(11)
φ
x
=
-
∂
w
0
∂
x
+
(
c
1
2
H
11
-
c
1
F
11
A
55
-
6
c
1
D
55
+
9
c
1
2
F
55
hhhhhhhhh
-
D
11
-
2
c
1
F
11
+
c
1
2
H
11
A
55
-
6
c
1
D
55
+
9
c
1
2
F
55
)
∂
3
w
0
∂
x
3
.In order to obtain the dimensionless equations, the variables and parameters are introduced as given below(12)
w
0
*
=
w
0
h
,
x
*
=
x
L
,
t
*
=
t
Q
11
I
I
0
L
4
,
A
11
=
Q
11
h
A
11
*
,
D
11
=
Q
11
h
3
3
D
11
*
,
D
55
=
Q
11
h
3
3
D
55
*
,
F
11
=
Q
11
h
5
5
F
11
*
,
H
11
=
Q
11
h
7
7
H
11
*
,
A
55
=
Q
11
h
A
55
*
,
f
*
=
L
4
h
Q
11
f
,
F
55
=
Q
11
h
5
5
F
55
*
.Naturally, for simplicity, the star is dropped in the following analysis.Closed form solutions forw
0 can be obtained for the beam by assuming
(13)
w
0
(
x
,
t
)
=
∑
w
k
(
t
)
sin
k
π
x
.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 inw
0 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
(14)
w
¨
1
+
c
3
w
˙
1
+
(
ω
2
-
π
2
f
)
w
1
+
3
π
4
D
11
8
w
1
3
=
0
,
where
(15)
ω
2
=
D
11
π
4
4
+
H
11
π
4
4
-
5
F
11
π
4
12
-
π
6
h
H
11
k
0
48
L
2
,
k
0
=
-
3
h
2
D
11
h
4
A
55
-
18
c
1
h
2
D
55
+
45
c
1
2
F
55
.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,
(16)
f
(
t
)
=
f
^
0
cos
⁡
(
Ω
t
+
γ
W
(
t
)
+
φ
)
,
where f
^
0 and Ω are the amplitude and center frequency of the excitation, which are constants, γ
≥
0 is an intensity, which represents the bandwidth, of the narrow-band random excitation, W
(
t
) is a standard Wiener process, and φ is a uniformly distributed random number in [
0,2
π
]. The inclusion of the phase angle φ makes f
(
t
) a stationary process. So f
(
t
) is a non-Gaussian distributed stationary stochastic process and has density function p
(
f
)
=
1
/
(
π
f
^
0
2
-
f
2
), zero mean value, and spectral density function as follows:
(17)
S
f
(
ω
)
=
1
4
π
f
^
0
2
γ
2
(
Ω
2
+
ω
2
+
γ
2
/
4
)
(
Ω
2
-
ω
2
+
γ
4
/
4
)
2
+
ω
2
γ
4
.The bandwidth of processf
(
t
) 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
(18)
w
¨
1
+
2
ς
^
ω
w
˙
1
+
ω
2
[
1
-
f
0
cos
⁡
(
Ω
t
+
γ
W
(
t
)
+
φ
)
]
w
1
+
α
^
w
1
3
=
0
,
where c
3
=
2
ς
^
ω, f
0
=
π
2
f
^
0
/
ω
2, α
^
=
3
π
4
D
11
/
8, and 0
<
ε
≪
1 is a small parameter.Also, a new time scale is defined asτ
=
ω
t. Let B
=
Ω
/
ω, ε
ς
=
ς
^, ε
α
=
α
^
/
ω
2, and ε
f
~
=
f
0 and rewrite (18) as
(19)
w
1
τ
τ
+
2
ε
ς
w
1
τ
+
[
1
-
ε
f
~
cos
⁡
(
B
τ
+
γ
W
(
τ
)
+
φ
)
]
w
1
+
ε
α
w
1
3
=
0
.


### 2.2. Analysis of Principal Parametric Resonance

In what follows, calculations have been performed for a graphite/epoxy symmetric laminated cross-ply beam with dimensionsL
=
0.5 m, b
=
0.1 m, h
=
0.05 m, ρ
=
1000 kg/m3, E
1
=
40 GPa, E
2
=
1 GPa, G
12
=
G
13
=
0.6 GPa, G
23
=
0.5 GPa, ν
12
=
0.25, and N
=
4
(
0
°
/
9
0
°
/
9
0
°
/
0
°
). Thus, parameter ε
α in (19) is calculated to be 2.4566 when the fiber orientation angle θ
1 of the first lamina with respect to the x-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(20)
w
1
(
T
0
,
T
1
;
ε
)
=
w
11
(
T
0
,
T
1
,
…
)
+
ε
w
12
(
T
0
,
T
1
,
…
)
+
⋯
.Moreover, for unit Wiener progressW
(
τ
), because of E
[
W
(
τ
)
]
=
0 and E
[
W
2
(
τ
)
]
=
τ, we have
(21)
γ
W
(
τ
)
=
(
γ
ε
)
W
(
ε
τ
)
=
γ
-
W
(
T
1
)
.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
(22)
B
=
2
+
ε
σ
.
Substituting (20) and (22) into (19) yieldsorderε
0:
(23)
D
0
2
w
11
+
w
11
=
0
,
orderε:
(24)
D
0
2
w
12
+
w
12
=
-
2
D
0
D
1
w
11
-
2
ς
D
0
w
11
+
f
~
w
11
cos
⁡
[
B
T
0
+
γ
-
W
(
T
1
)
+
φ
]
-
α
w
11
3
.The first-order approximate solution of (23) is
(25)
w
11
(
T
0
,
T
1
)
=
A
(
T
1
)
e
i
T
0
+
cc
,
where cc stands for the complex conjugate of the preceding terms.Substituting (22) and (25) into (24) and eliminating the secular terms, we have
(26)
2
i
(
D
1
A
+
ς
A
)
-
f
~
2
A
-
e
i
(
σ
T
1
+
γ
-
W
+
φ
)
+
α
e
A
2
A
-
=
0
,
where α
e
=
3
α.Substituting the polar form(27)
A
(
T
1
)
=
1
2
a
(
T
1
)
e
i
(
σ
T
1
+
γ
-
W
+
φ
-
β
)
/
2
,
into (26) and separating the real and imaginary parts, we get(28a)
D
1
a
=
-
ς
a
+
f
~
4
a
sin
β
,
(28b)
a
D
1
β
=
σ
a
+
f
~
2
a
cos
⁡
β
-
1
4
α
e
a
3
+
a
γ
-
D
1
W
.


#### 2.2.1. Stability Analysis of the Trivial Response

It can be found that (28a) and (28b) have a trivial steady state solution a
=
0. In order to determine the stability of the trivial steady state response, the trivial response is expanded near a
=
0 and the corresponding nonlinear terms are neglected; the linearization form of (28a) and (28b) can finally be expressed as(29a)
D
1
a
=
-
ς
a
+
f
~
4
a
sin
β
,
(29b)
D
1
β
=
σ
+
f
~
2
cos
⁡
β
+
γ
-
D
1
W
.Letr
=
ln
⁡
a and we rewrite (29a) and (29b) in the form of Itô equation as given below(30a)
d
r
=
(
-
ς
+
f
~
4
sin
β
)
d
T
1
,
(30b)
d
β
=
(
σ
+
f
~
2
cos
⁡
β
)
d
T
1
+
γ
-
d
W
.It can be found that the stochastic processβ generated on [
0,2
π
] by (30a) and (30b) is Markov and is also ergodic on [
0,2
π
] since the diffusion process is nonsingular. Naturally, the invariant measure, that is, the steady state probability density function p
(
β
) of the process β, can be determined by the following FPK equation:
(31)
d
2
p
(
β
)
d
β
2
-
d
d
β
[
(
σ
-
+
f
-
cos
⁡
β
)
p
(
β
)
]
=
0
,
where σ
-
=
2
σ
/
γ
-
2 and f
-
=
f
~
/
γ
-
2.The unique solution of (31) satisfies both the periodicity condition p
(
β
)
=
p
(
β
+
2
π
) and the normality condition ∫
0
2
π
p
(
β
)
d
β
=
1, respectively. Thus, the solution of (31) can be expressed as
(32)
p
(
β
)
=
1
C
e
(
σ
-
β
+
f
-
sin
β
)
∫
β
β
+
2
π
e
[
-
(
σ
-
x
+
f
-
sin
x
)
]
d
x
,
where C
=
4
π
2
e
(
-
σ
-
π
)
I
i
σ
-
(
-
f
-
)
I
-
i
σ
-
(
-
f
-
) and I
n
(
x
) 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
(33)
λ
(
a
0
,
β
0
)
=
lim
⁡
T
1
→
∞
1
T
1
ln
⁡
|
a
(
T
;
a
0
,
β
0
)
|
w
.
p
.
1
,
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 λ
=
λ
max
⁡; that is, λ
=
0 is the bifurcation point of the stability of the trivial response. According to [30], the largest Lyapunov exponent λ is given as
(34)
λ
=
lim
⁡
T
1
→
∞
1
T
1
ln
⁡
|
a
(
T
1
)
a
(
0
)
|
=
-
ς
+
lim
⁡
T
1
→
∞
1
T
1
∫
0
T
1
f
~
4
sin
β
(
τ
)
d
τ
=
-
ς
+
f
~
4
E
[
sin
β
]
=
-
ς
+
f
~
4
∫
0
2
π
sin
β
p
(
β
)
d
β
=
-
ς
+
f
~
4
[
I
1
-
i
σ
-
(
f
-
)
I
-
i
σ
-
(
f
-
)
+
I
1
+
i
σ
-
(
f
-
)
I
i
σ
-
(
f
-
)
]
.The variation of the largest Lyapunov exponentλ as λ
-
(
B
,
f
0
) 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.01, 0.1, and 1 and ε
ς
=
0.1, 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 B
=
2, 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 λ
=
0 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 f
0 but intensify the instability of the trivial response for a higher amplitude f
0. Such results can also be found in Figure 5 for three isohypse curves of λ
=
0 resulting from different γ; that is, γ
=
0.01, γ
=
0.1, and γ
=
1.0, respectively.The largest Lyapunov exponentλ of the trivial response of the system: γ
=
0.01, ε
ς
=
0.1. (a) Mesh surface of λ
-
(
B
,
f
0
); (b) isohypse curves of λ.
(a)
(b)The largest Lyapunov exponentλ of the trivial response of the system: γ
=
0.1, ε
ς
=
0.1. (a) Mesh surface of λ
-
(
B
,
f
0
); (b) isohypse curves of λ.
(a)
(b)The largest Lyapunov exponentλ of the trivial response of the system: γ
=
1.0, ε
ς
=
0.1. (a) Mesh surface of λ
-
(
B
,
f
0
); (b) isohypse curves of λ.
(a)
(b)Figure 5
Isohypse curves ofλ
=
0.


#### 2.2.2. Stochastic Jump and Bifurcation

As the reduced (28b) has the coupled term a
D
1
β, 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]
(35)
x
=
a
sin
β
,
y
=
a
cos
⁡
β
,
and, introducing it into (28a) and (28b), we have(36a)
D
1
x
=
-
ς
x
-
f
~
4
x
y
x
2
+
y
2
-
σ
y
+
1
4
α
e
y
(
x
2
+
y
2
)
-
y
γ
-
D
1
W
,
(36b)
D
1
y
=
-
ς
y
+
f
~
4
(
y
2
+
2
x
2
)
x
2
+
y
2
+
σ
x
-
1
4
α
e
x
(
x
2
+
y
2
)
+
x
γ
-
D
1
W
,where x and y are the two perpendicular components of a; that is, a
2
=
x
2
+
y
2.Equations (36a) and (36b) are equivalent to the following set of Itô equation:(37a)
d
x
=
m
1
(
x
,
y
,
T
1
)
d
T
1
+
σ
11
(
x
,
y
,
T
1
)
d
W
+
σ
12
(
x
,
y
,
T
1
)
d
W
,
(37b)
d
y
=
m
2
(
x
,
y
,
T
1
)
d
T
1
+
σ
21
(
x
,
y
,
T
1
)
d
W
+
σ
22
(
x
,
y
,
T
1
)
d
W
,where
(38)
m
1
(
x
,
y
,
T
1
)
=
-
ς
x
-
f
~
4
x
y
x
2
+
y
2
-
σ
y
+
1
4
α
e
y
(
x
2
+
y
2
)
-
1
2
γ
-
2
x
,
b
11
(
x
,
y
,
T
1
)
=
y
2
γ
-
2
,
b
22
(
x
,
y
,
T
1
)
=
x
2
γ
-
2
,
m
2
(
x
,
y
,
T
1
)
=
-
ς
y
+
f
~
4
(
y
2
+
2
x
2
)
x
2
+
y
2
+
σ
x
-
1
4
α
e
x
(
x
2
+
y
2
)
-
1
2
γ
-
2
y
,
b
12
(
x
,
y
,
T
1
)
=
b
21
(
x
,
y
,
T
1
)
=
-
x
y
γ
-
2
,
b
i
j
(
x
,
y
,
T
1
)
=
∑
l
=
1
2
σ
i
l
(
X
,
Y
,
T
1
)
σ
j
l
(
X
,
Y
,
T
1
)
|
X
=
x
Y
=
y
,
h
=
h
h
h
h
h
h
h
h
h
h
h
h
h
h
i
,
j
=
1,2
,
and the last terms in m
1 and m
2 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
(39)
∂
p
∂
T
1
=
-
∂
∂
x
[
m
1
p
]
-
∂
∂
y
[
m
2
p
]
+
1
2
∂
2
(
b
11
p
)
∂
x
2
+
1
2
∂
2
(
b
22
p
)
∂
y
2
+
1
2
∂
2
(
b
12
p
)
∂
x
∂
y
+
1
2
∂
2
(
b
21
p
)
∂
y
∂
x
,
where p
=
p
(
x
,
y
,
T
1
∣
x
0
,
y
0
) is the transition probability density of x and y.The initial condition of FPK equation (39) is
(40)
p
=
δ
(
x
-
x
0
)
δ
(
y
-
y
0
)
,
T
1
=
0
.
The boundary conditions with respect to x and y are
(41)
p
,
∂
p
∂
x
,
∂
p
∂
y
⟶
0
,
as
|
x
|
,
|
y
|
⟶
∞
.Whenγ
-
=
0 and φ
=
0, (28a) and (28b) is a deterministically excited system and it has one trivial solution and two possible steady state nontrivial solutions when f
~
>
4
ς as given below
(42)
a
0
=
4
(
B
-
2
)
3
ε
α
±
2
(
ε
f
~
)
2
-
(
4
ε
ς
)
2
3
ε
α
.The theoretical frequency-response curves based on (42) are shown in Figure 6 for the aforementioned beam when γ
-
=
0 and φ
=
0. 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 whenf
0
=
1 and ε
ς
=
0.1.Figure7 lists a series of change of the stationary joint probability density of x and y, which is numerically solved from the FPK equation (39) by using finite difference method [34, 35] to different excitation central frequency B when γ
=
0.12 and f
0
=
1. Numerical investigation shows that the numerical result becomes convergent when the grid number is greater than 65
×
65 within the amplitude region [
(
-
2
,
2
)
,
(
-
2
,
2
)
]; thus we choose the grid number M
×
N
=
71
×
71 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 B 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 B 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 B 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 frequenciesB: γ
=
0.12, f
0
=
1. (a) B
=
2.3; (b) B
=
3.0; (c) B
=
3.5; (d) B
=
4.0; (e) B
=
5.0; (f) B
=
6.0.
(a)
B
=
2.3
(b)
B
=
3.0
(c)
B
=
3.5
(d)
B
=
4.0
(e)
B
=
5.0
(f)
B
=
6.0Figure8 reveals the stationary joint probability density to different bandwidths when B
=
2.8 and f
0
=
1.0. Hear the choice of B
=
2.8 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, γ
=
0.12 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)).Stationary joint probability density of the system to different bandwidthsγ. B
=
2.8, f
0
=
1.
(a)
γ
=
0.12
(b)
γ
=
0.20
(c)
γ
=
0.30
(d)
γ
=
0.50The 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 whenB
=
2.7.Figure10 shows the stationary joint probability density of the system to different excitation forces when B
=
2.7 and γ
=
0.10. We can find that there is only an outer flabellate peak when the excitation force, for instance, f
0
=
2.0, 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.B
=
2.7, γ
=
0.10.
(a)
f
0
=
2.0
(b)
f
0
=
1.4
(c)
f
0
=
1.0
(d)
f
0
=
0.6In order to bring forth the influence of the layup or laminate properties, in what follows, the fiber orientation angleθ
1 of the first lamina with respect to the x-axis of the beam is taken into consideration. Table 1 lists a series of data for the dimensionless natural frequency ω, the dimensionless excitation amplitude f
0, the dimensionless central frequency B, and the nonlinear parameter ε
α with a smaller change of the fiber orientation angle θ
1. 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 f
0
ω
2, and the bandwidth γ in (18) are fixed; here we choose Ω
=
5.4200, f
0
ω
2
=
3.2642, and γ
=
0.12.Table 1
System and excitation parameters corresponding toθ
1.
θ
1
ω
α
^
ε
α
f
0
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.2324Figure11 shows the stationary joint probability density of the system to different fiber orientation angles when Ω
=
5.4200, f
0
ω
2
=
3.2642, and γ
=
0.12. 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 θ
1
≤
2
° (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 θ
1
=
3
°, the probable motion slightly moves towards the stable stationary nontrivial branch (see Figure 11(d)). As the fiber orientation angle increases from 4
° to 9
°, 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 (1
0
°
-
1
1
°), the probable motion goes back towards the stable stationary nontrivial branch (see Figures 11(k)-11(l)). When θ
1
=
1
2
°, 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 θ
1
=
1
6
° (see Figure 11(q)).Stationary joint probability density of the system to different fiber orientation anglesθ
1 of the first lamina with respect to the x-axis of the beam. Ω
=
5.4200, f
0
ω
2
=
3.2642, and γ
=
0.12.
(a)
θ
1
=
0
°
(b)
θ
1
=
1
°
(c)
θ
1
=
2
°
(d)
θ
1
=
3
°
(e)
θ
1
=
4
°
(f)
θ
1
=
5
°
(g)
θ
1
=
6
°
(h)
θ
1
=
7
°
(i)
θ
1
=
8
°
(j)
θ
1
=
9
°
(k)
θ
1
=
1
0
°
(l)
θ
1
=
1
1
°
(m)
θ
1
=
1
2
°
(n)
θ
1
=
1
3
°
(o)
θ
1
=
1
4
°
(p)
θ
1
=
1
5
°
(q)
θ
1
=
1
6
°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. Figures12(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 x and y in both Figures 11(f) and 11(p) to the amplitude a, 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θ
1 of the first lamina with respect to the x-axis of the beam, Ω
=
5.4200, f
0
ω
2
=
3.2642, and γ
=
0.12.
(a)
θ
1
=
5
°
(b)
θ
1
=
1
5
°


## 2.1. Equation of Motion

A simply supported laminated composite beam is shown in Figure1 with length L and thickness h subject to axial load f. The displacement proposed by Levinson [4] and used by Reddy [6] is given as below(1a)
u
(
x
,
z
,
t
)
=
u
0
(
x
,
t
)
+
z
φ
x
(
x
,
t
)
-
z
3
c
1
(
φ
x
+
∂
w
0
∂
x
)
,
(1b)
w
(
x
,
z
,
t
)
=
w
0
(
x
,
t
)
,
(1c)
φ
x
=
(
∂
u
∂
z
)
z
=
0
,where c
1
=
4
/
3
h
2, u and w are total displacements, u
0 and w
0 are the displacements of a point in the middle surface of the beam in the x and z directions, and φ
x is rotation of the normal about the y-axis.Figure 1
Configuration of a laminated composite beam.According to Reddy and von Karman’s theory, the governing equations of the motion by the Hamilton principle can be obtained as(2a)
A
11
∂
2
u
0
∂
x
2
+
A
11
∂
w
0
∂
x
∂
2
w
0
∂
x
2
+
(
B
11
-
c
1
E
11
)
∂
2
φ
x
∂
x
2
-
c
1
E
11
∂
3
w
0
∂
x
3
=
I
0
∂
2
u
0
∂
t
2
+
J
1
∂
2
φ
x
∂
t
2
-
c
1
I
3
∂
3
w
0
∂
x
∂
t
2
,
(2b)
A
11
∂
u
0
∂
x
∂
2
w
0
∂
x
2
+
A
11
∂
2
u
0
∂
x
2
∂
w
0
∂
x
+
c
1
E
11
∂
3
u
0
∂
x
3
+
(
A
55
-
6
c
1
D
55
+
9
c
1
2
F
55
)
∂
φ
x
∂
x
+
3
2
A
11
(
∂
w
0
∂
x
)
2
∂
2
w
0
∂
x
2
+
(
B
11
-
c
1
E
11
)
∂
2
φ
x
∂
x
2
∂
w
0
∂
x
+
c
1
(
F
11
-
c
1
H
11
)
∂
3
φ
x
∂
x
3
+
(
A
55
-
6
c
1
D
55
+
9
c
1
2
F
55
)
∂
2
w
0
∂
x
2
+
(
B
11
-
c
1
E
11
)
×
∂
φ
x
∂
x
∂
2
w
0
∂
x
2
-
c
1
2
H
11
∂
4
w
0
∂
x
4
-
c
3
w
˙
0
=
I
0
∂
2
w
0
∂
t
2
+
c
1
I
3
∂
3
u
0
∂
x
∂
t
2
+
c
1
J
4
∂
3
φ
x
∂
x
∂
t
2
-
c
1
2
I
6
∂
4
w
0
∂
x
2
∂
t
2
,
(2c)
(
B
11
-
c
1
E
11
)
∂
2
u
0
∂
x
2
+
(
B
11
-
c
1
E
11
)
∂
w
0
∂
x
∂
2
w
0
∂
x
2
-
(
c
1
F
11
-
c
1
2
H
11
)
∂
3
w
0
∂
x
3
+
[
(
D
11
-
c
1
F
11
)
-
c
1
(
F
11
-
c
1
H
11
)
]
∂
2
φ
x
∂
x
2
-
(
A
55
-
6
c
1
D
55
+
9
c
1
2
F
55
)
φ
x
-
(
A
55
-
6
c
1
D
55
+
9
c
1
2
F
55
)
∂
w
0
∂
x
=
J
1
∂
2
u
0
∂
t
2
+
K
2
∂
2
φ
x
∂
t
2
-
c
1
J
4
∂
3
w
0
∂
x
∂
t
2
,where c
3 is the coefficient of damping and other coefficients can be found in the appendix.Also, the boundary conditions can be written asx
=
0:(3a)
δ
u
0
=
0
or
b
(
A
11
∂
u
0
∂
x
+
1
2
A
11
(
∂
w
0
∂
x
)
2
h
h
h
(
∂
w
0
∂
x
)
2
+
(
B
11
-
c
1
E
11
)
∂
φ
x
∂
x
-
c
1
E
11
∂
2
w
0
∂
x
2
)
+
f
h
=
0
,
x
=
L:
(3b)
δ
u
0
=
0
or
A
11
∂
u
0
∂
x
+
1
2
A
11
(
∂
w
0
∂
x
)
2
+
(
B
11
-
c
1
E
11
)
∂
φ
x
∂
x
-
c
1
E
11
∂
2
w
0
∂
x
2
=
0
,
x
=
0 and x
=
L:(4)
δ
φ
x
=
0
or
(
B
11
z
-
c
1
E
11
z
3
)
[
∂
u
0
∂
x
+
1
2
(
∂
w
0
∂
x
)
2
]
+
(
D
11
-
2
c
1
F
11
+
c
1
2
H
11
)
∂
φ
x
∂
x
-
(
c
1
F
11
-
c
1
2
H
11
)
∂
2
w
0
∂
x
2
=
0
,
x
=
0 and x
=
L:
(5)
∂
δ
w
0
∂
x
=
0
or
E
11
∂
u
0
∂
x
+
1
2
E
11
(
∂
w
0
∂
x
)
2
+
c
1
(
F
11
-
c
1
H
11
)
∂
φ
x
∂
x
-
c
1
2
H
11
∂
2
w
0
∂
x
2
=
0
,
x
=
0 and x
=
L:
(6)
δ
w
0
=
0
,
or
[
A
11
∂
u
0
∂
x
+
1
2
A
11
(
∂
w
0
∂
x
)
2
h
.
+
(
B
11
-
c
1
E
11
)
∂
φ
x
∂
x
-
c
1
E
11
∂
2
w
0
∂
x
2
]
∂
w
0
∂
x
+
c
1
[
E
11
∂
3
u
0
∂
x
3
+
E
11
(
∂
2
w
0
∂
x
2
)
2
+
E
11
∂
w
0
∂
x
∂
3
w
0
∂
x
3
h
.
h
.
h
.
h
.
+
(
F
11
-
c
1
H
11
)
∂
3
φ
x
∂
x
3
-
c
1
H
11
∂
4
w
0
∂
x
4
]
+
(
A
55
-
6
c
1
D
55
+
9
c
1
2
F
55
)
(
φ
x
+
∂
w
0
∂
x
)
+
c
1
I
3
∂
2
u
0
∂
t
2
-
c
1
(
I
4
-
c
1
I
6
)
∂
2
φ
x
∂
t
2
+
c
1
2
I
6
∂
3
w
0
∂
x
∂
t
2
=
0
.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(7)
B
11
=
E
11
=
0
,
I
1
=
I
3
=
0
.Substituting (7) into (2a) and (2c) and neglecting the longitudinal inertia and the moment of inertia about the neutral axis, we have
(8)
A
11
∂
2
u
0
∂
x
2
+
A
11
∂
w
0
∂
x
∂
2
w
0
∂
x
2
=
0
,
(9)
-
(
A
55
-
6
c
1
D
55
+
9
c
1
2
F
55
)
∂
w
0
∂
x
-
(
A
55
-
6
c
1
D
55
+
9
c
1
2
F
55
)
φ
x
+
[
(
D
11
-
c
1
F
11
)
-
c
1
(
F
11
-
c
1
H
11
)
]
∂
2
φ
x
∂
x
2
-
(
c
1
F
11
-
c
1
2
H
11
)
∂
3
w
0
∂
x
3
=
0
.Integrating (8) and using (3a) and (3b), we arrive at
(10)
∂
u
0
∂
x
=
-
1
2
(
∂
w
0
∂
x
)
2
+
h
f
A
11
.According to the boundary conditions, the solution of (9) can be approximately expressed as [33]
(11)
φ
x
=
-
∂
w
0
∂
x
+
(
c
1
2
H
11
-
c
1
F
11
A
55
-
6
c
1
D
55
+
9
c
1
2
F
55
hhhhhhhhh
-
D
11
-
2
c
1
F
11
+
c
1
2
H
11
A
55
-
6
c
1
D
55
+
9
c
1
2
F
55
)
∂
3
w
0
∂
x
3
.In order to obtain the dimensionless equations, the variables and parameters are introduced as given below(12)
w
0
*
=
w
0
h
,
x
*
=
x
L
,
t
*
=
t
Q
11
I
I
0
L
4
,
A
11
=
Q
11
h
A
11
*
,
D
11
=
Q
11
h
3
3
D
11
*
,
D
55
=
Q
11
h
3
3
D
55
*
,
F
11
=
Q
11
h
5
5
F
11
*
,
H
11
=
Q
11
h
7
7
H
11
*
,
A
55
=
Q
11
h
A
55
*
,
f
*
=
L
4
h
Q
11
f
,
F
55
=
Q
11
h
5
5
F
55
*
.Naturally, for simplicity, the star is dropped in the following analysis.Closed form solutions forw
0 can be obtained for the beam by assuming
(13)
w
0
(
x
,
t
)
=
∑
w
k
(
t
)
sin
k
π
x
.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 inw
0 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
(14)
w
¨
1
+
c
3
w
˙
1
+
(
ω
2
-
π
2
f
)
w
1
+
3
π
4
D
11
8
w
1
3
=
0
,
where
(15)
ω
2
=
D
11
π
4
4
+
H
11
π
4
4
-
5
F
11
π
4
12
-
π
6
h
H
11
k
0
48
L
2
,
k
0
=
-
3
h
2
D
11
h
4
A
55
-
18
c
1
h
2
D
55
+
45
c
1
2
F
55
.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,
(16)
f
(
t
)
=
f
^
0
cos
⁡
(
Ω
t
+
γ
W
(
t
)
+
φ
)
,
where f
^
0 and Ω are the amplitude and center frequency of the excitation, which are constants, γ
≥
0 is an intensity, which represents the bandwidth, of the narrow-band random excitation, W
(
t
) is a standard Wiener process, and φ is a uniformly distributed random number in [
0,2
π
]. The inclusion of the phase angle φ makes f
(
t
) a stationary process. So f
(
t
) is a non-Gaussian distributed stationary stochastic process and has density function p
(
f
)
=
1
/
(
π
f
^
0
2
-
f
2
), zero mean value, and spectral density function as follows:
(17)
S
f
(
ω
)
=
1
4
π
f
^
0
2
γ
2
(
Ω
2
+
ω
2
+
γ
2
/
4
)
(
Ω
2
-
ω
2
+
γ
4
/
4
)
2
+
ω
2
γ
4
.The bandwidth of processf
(
t
) 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
(18)
w
¨
1
+
2
ς
^
ω
w
˙
1
+
ω
2
[
1
-
f
0
cos
⁡
(
Ω
t
+
γ
W
(
t
)
+
φ
)
]
w
1
+
α
^
w
1
3
=
0
,
where c
3
=
2
ς
^
ω, f
0
=
π
2
f
^
0
/
ω
2, α
^
=
3
π
4
D
11
/
8, and 0
<
ε
≪
1 is a small parameter.Also, a new time scale is defined asτ
=
ω
t. Let B
=
Ω
/
ω, ε
ς
=
ς
^, ε
α
=
α
^
/
ω
2, and ε
f
~
=
f
0 and rewrite (18) as
(19)
w
1
τ
τ
+
2
ε
ς
w
1
τ
+
[
1
-
ε
f
~
cos
⁡
(
B
τ
+
γ
W
(
τ
)
+
φ
)
]
w
1
+
ε
α
w
1
3
=
0
.


## 2.2. Analysis of Principal Parametric Resonance

In what follows, calculations have been performed for a graphite/epoxy symmetric laminated cross-ply beam with dimensionsL
=
0.5 m, b
=
0.1 m, h
=
0.05 m, ρ
=
1000 kg/m3, E
1
=
40 GPa, E
2
=
1 GPa, G
12
=
G
13
=
0.6 GPa, G
23
=
0.5 GPa, ν
12
=
0.25, and N
=
4
(
0
°
/
9
0
°
/
9
0
°
/
0
°
). Thus, parameter ε
α in (19) is calculated to be 2.4566 when the fiber orientation angle θ
1 of the first lamina with respect to the x-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(20)
w
1
(
T
0
,
T
1
;
ε
)
=
w
11
(
T
0
,
T
1
,
…
)
+
ε
w
12
(
T
0
,
T
1
,
…
)
+
⋯
.Moreover, for unit Wiener progressW
(
τ
), because of E
[
W
(
τ
)
]
=
0 and E
[
W
2
(
τ
)
]
=
τ, we have
(21)
γ
W
(
τ
)
=
(
γ
ε
)
W
(
ε
τ
)
=
γ
-
W
(
T
1
)
.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
(22)
B
=
2
+
ε
σ
.
Substituting (20) and (22) into (19) yieldsorderε
0:
(23)
D
0
2
w
11
+
w
11
=
0
,
orderε:
(24)
D
0
2
w
12
+
w
12
=
-
2
D
0
D
1
w
11
-
2
ς
D
0
w
11
+
f
~
w
11
cos
⁡
[
B
T
0
+
γ
-
W
(
T
1
)
+
φ
]
-
α
w
11
3
.The first-order approximate solution of (23) is
(25)
w
11
(
T
0
,
T
1
)
=
A
(
T
1
)
e
i
T
0
+
cc
,
where cc stands for the complex conjugate of the preceding terms.Substituting (22) and (25) into (24) and eliminating the secular terms, we have
(26)
2
i
(
D
1
A
+
ς
A
)
-
f
~
2
A
-
e
i
(
σ
T
1
+
γ
-
W
+
φ
)
+
α
e
A
2
A
-
=
0
,
where α
e
=
3
α.Substituting the polar form(27)
A
(
T
1
)
=
1
2
a
(
T
1
)
e
i
(
σ
T
1
+
γ
-
W
+
φ
-
β
)
/
2
,
into (26) and separating the real and imaginary parts, we get(28a)
D
1
a
=
-
ς
a
+
f
~
4
a
sin
β
,
(28b)
a
D
1
β
=
σ
a
+
f
~
2
a
cos
⁡
β
-
1
4
α
e
a
3
+
a
γ
-
D
1
W
.


### 2.2.1. Stability Analysis of the Trivial Response

It can be found that (28a) and (28b) have a trivial steady state solution a
=
0. In order to determine the stability of the trivial steady state response, the trivial response is expanded near a
=
0 and the corresponding nonlinear terms are neglected; the linearization form of (28a) and (28b) can finally be expressed as(29a)
D
1
a
=
-
ς
a
+
f
~
4
a
sin
β
,
(29b)
D
1
β
=
σ
+
f
~
2
cos
⁡
β
+
γ
-
D
1
W
.Letr
=
ln
⁡
a and we rewrite (29a) and (29b) in the form of Itô equation as given below(30a)
d
r
=
(
-
ς
+
f
~
4
sin
β
)
d
T
1
,
(30b)
d
β
=
(
σ
+
f
~
2
cos
⁡
β
)
d
T
1
+
γ
-
d
W
.It can be found that the stochastic processβ generated on [
0,2
π
] by (30a) and (30b) is Markov and is also ergodic on [
0,2
π
] since the diffusion process is nonsingular. Naturally, the invariant measure, that is, the steady state probability density function p
(
β
) of the process β, can be determined by the following FPK equation:
(31)
d
2
p
(
β
)
d
β
2
-
d
d
β
[
(
σ
-
+
f
-
cos
⁡
β
)
p
(
β
)
]
=
0
,
where σ
-
=
2
σ
/
γ
-
2 and f
-
=
f
~
/
γ
-
2.The unique solution of (31) satisfies both the periodicity condition p
(
β
)
=
p
(
β
+
2
π
) and the normality condition ∫
0
2
π
p
(
β
)
d
β
=
1, respectively. Thus, the solution of (31) can be expressed as
(32)
p
(
β
)
=
1
C
e
(
σ
-
β
+
f
-
sin
β
)
∫
β
β
+
2
π
e
[
-
(
σ
-
x
+
f
-
sin
x
)
]
d
x
,
where C
=
4
π
2
e
(
-
σ
-
π
)
I
i
σ
-
(
-
f
-
)
I
-
i
σ
-
(
-
f
-
) and I
n
(
x
) 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
(33)
λ
(
a
0
,
β
0
)
=
lim
⁡
T
1
→
∞
1
T
1
ln
⁡
|
a
(
T
;
a
0
,
β
0
)
|
w
.
p
.
1
,
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 λ
=
λ
max
⁡; that is, λ
=
0 is the bifurcation point of the stability of the trivial response. According to [30], the largest Lyapunov exponent λ is given as
(34)
λ
=
lim
⁡
T
1
→
∞
1
T
1
ln
⁡
|
a
(
T
1
)
a
(
0
)
|
=
-
ς
+
lim
⁡
T
1
→
∞
1
T
1
∫
0
T
1
f
~
4
sin
β
(
τ
)
d
τ
=
-
ς
+
f
~
4
E
[
sin
β
]
=
-
ς
+
f
~
4
∫
0
2
π
sin
β
p
(
β
)
d
β
=
-
ς
+
f
~
4
[
I
1
-
i
σ
-
(
f
-
)
I
-
i
σ
-
(
f
-
)
+
I
1
+
i
σ
-
(
f
-
)
I
i
σ
-
(
f
-
)
]
.The variation of the largest Lyapunov exponentλ as λ
-
(
B
,
f
0
) 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.01, 0.1, and 1 and ε
ς
=
0.1, 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 B
=
2, 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 λ
=
0 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 f
0 but intensify the instability of the trivial response for a higher amplitude f
0. Such results can also be found in Figure 5 for three isohypse curves of λ
=
0 resulting from different γ; that is, γ
=
0.01, γ
=
0.1, and γ
=
1.0, respectively.The largest Lyapunov exponentλ of the trivial response of the system: γ
=
0.01, ε
ς
=
0.1. (a) Mesh surface of λ
-
(
B
,
f
0
); (b) isohypse curves of λ.
(a)
(b)The largest Lyapunov exponentλ of the trivial response of the system: γ
=
0.1, ε
ς
=
0.1. (a) Mesh surface of λ
-
(
B
,
f
0
); (b) isohypse curves of λ.
(a)
(b)The largest Lyapunov exponentλ of the trivial response of the system: γ
=
1.0, ε
ς
=
0.1. (a) Mesh surface of λ
-
(
B
,
f
0
); (b) isohypse curves of λ.
(a)
(b)Figure 5
Isohypse curves ofλ
=
0.


### 2.2.2. Stochastic Jump and Bifurcation

As the reduced (28b) has the coupled term a
D
1
β, 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]
(35)
x
=
a
sin
β
,
y
=
a
cos
⁡
β
,
and, introducing it into (28a) and (28b), we have(36a)
D
1
x
=
-
ς
x
-
f
~
4
x
y
x
2
+
y
2
-
σ
y
+
1
4
α
e
y
(
x
2
+
y
2
)
-
y
γ
-
D
1
W
,
(36b)
D
1
y
=
-
ς
y
+
f
~
4
(
y
2
+
2
x
2
)
x
2
+
y
2
+
σ
x
-
1
4
α
e
x
(
x
2
+
y
2
)
+
x
γ
-
D
1
W
,where x and y are the two perpendicular components of a; that is, a
2
=
x
2
+
y
2.Equations (36a) and (36b) are equivalent to the following set of Itô equation:(37a)
d
x
=
m
1
(
x
,
y
,
T
1
)
d
T
1
+
σ
11
(
x
,
y
,
T
1
)
d
W
+
σ
12
(
x
,
y
,
T
1
)
d
W
,
(37b)
d
y
=
m
2
(
x
,
y
,
T
1
)
d
T
1
+
σ
21
(
x
,
y
,
T
1
)
d
W
+
σ
22
(
x
,
y
,
T
1
)
d
W
,where
(38)
m
1
(
x
,
y
,
T
1
)
=
-
ς
x
-
f
~
4
x
y
x
2
+
y
2
-
σ
y
+
1
4
α
e
y
(
x
2
+
y
2
)
-
1
2
γ
-
2
x
,
b
11
(
x
,
y
,
T
1
)
=
y
2
γ
-
2
,
b
22
(
x
,
y
,
T
1
)
=
x
2
γ
-
2
,
m
2
(
x
,
y
,
T
1
)
=
-
ς
y
+
f
~
4
(
y
2
+
2
x
2
)
x
2
+
y
2
+
σ
x
-
1
4
α
e
x
(
x
2
+
y
2
)
-
1
2
γ
-
2
y
,
b
12
(
x
,
y
,
T
1
)
=
b
21
(
x
,
y
,
T
1
)
=
-
x
y
γ
-
2
,
b
i
j
(
x
,
y
,
T
1
)
=
∑
l
=
1
2
σ
i
l
(
X
,
Y
,
T
1
)
σ
j
l
(
X
,
Y
,
T
1
)
|
X
=
x
Y
=
y
,
h
=
h
h
h
h
h
h
h
h
h
h
h
h
h
h
i
,
j
=
1,2
,
and the last terms in m
1 and m
2 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
(39)
∂
p
∂
T
1
=
-
∂
∂
x
[
m
1
p
]
-
∂
∂
y
[
m
2
p
]
+
1
2
∂
2
(
b
11
p
)
∂
x
2
+
1
2
∂
2
(
b
22
p
)
∂
y
2
+
1
2
∂
2
(
b
12
p
)
∂
x
∂
y
+
1
2
∂
2
(
b
21
p
)
∂
y
∂
x
,
where p
=
p
(
x
,
y
,
T
1
∣
x
0
,
y
0
) is the transition probability density of x and y.The initial condition of FPK equation (39) is
(40)
p
=
δ
(
x
-
x
0
)
δ
(
y
-
y
0
)
,
T
1
=
0
.
The boundary conditions with respect to x and y are
(41)
p
,
∂
p
∂
x
,
∂
p
∂
y
⟶
0
,
as
|
x
|
,
|
y
|
⟶
∞
.Whenγ
-
=
0 and φ
=
0, (28a) and (28b) is a deterministically excited system and it has one trivial solution and two possible steady state nontrivial solutions when f
~
>
4
ς as given below
(42)
a
0
=
4
(
B
-
2
)
3
ε
α
±
2
(
ε
f
~
)
2
-
(
4
ε
ς
)
2
3
ε
α
.The theoretical frequency-response curves based on (42) are shown in Figure 6 for the aforementioned beam when γ
-
=
0 and φ
=
0. 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 whenf
0
=
1 and ε
ς
=
0.1.Figure7 lists a series of change of the stationary joint probability density of x and y, which is numerically solved from the FPK equation (39) by using finite difference method [34, 35] to different excitation central frequency B when γ
=
0.12 and f
0
=
1. Numerical investigation shows that the numerical result becomes convergent when the grid number is greater than 65
×
65 within the amplitude region [
(
-
2
,
2
)
,
(
-
2
,
2
)
]; thus we choose the grid number M
×
N
=
71
×
71 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 B 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 B 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 B 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 frequenciesB: γ
=
0.12, f
0
=
1. (a) B
=
2.3; (b) B
=
3.0; (c) B
=
3.5; (d) B
=
4.0; (e) B
=
5.0; (f) B
=
6.0.
(a)
B
=
2.3
(b)
B
=
3.0
(c)
B
=
3.5
(d)
B
=
4.0
(e)
B
=
5.0
(f)
B
=
6.0Figure8 reveals the stationary joint probability density to different bandwidths when B
=
2.8 and f
0
=
1.0. Hear the choice of B
=
2.8 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, γ
=
0.12 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)).Stationary joint probability density of the system to different bandwidthsγ. B
=
2.8, f
0
=
1.
(a)
γ
=
0.12
(b)
γ
=
0.20
(c)
γ
=
0.30
(d)
γ
=
0.50The 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 whenB
=
2.7.Figure10 shows the stationary joint probability density of the system to different excitation forces when B
=
2.7 and γ
=
0.10. We can find that there is only an outer flabellate peak when the excitation force, for instance, f
0
=
2.0, 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.B
=
2.7, γ
=
0.10.
(a)
f
0
=
2.0
(b)
f
0
=
1.4
(c)
f
0
=
1.0
(d)
f
0
=
0.6In order to bring forth the influence of the layup or laminate properties, in what follows, the fiber orientation angleθ
1 of the first lamina with respect to the x-axis of the beam is taken into consideration. Table 1 lists a series of data for the dimensionless natural frequency ω, the dimensionless excitation amplitude f
0, the dimensionless central frequency B, and the nonlinear parameter ε
α with a smaller change of the fiber orientation angle θ
1. 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 f
0
ω
2, and the bandwidth γ in (18) are fixed; here we choose Ω
=
5.4200, f
0
ω
2
=
3.2642, and γ
=
0.12.Table 1
System and excitation parameters corresponding toθ
1.
θ
1
ω
α
^
ε
α
f
0
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.2324Figure11 shows the stationary joint probability density of the system to different fiber orientation angles when Ω
=
5.4200, f
0
ω
2
=
3.2642, and γ
=
0.12. 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 θ
1
≤
2
° (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 θ
1
=
3
°, the probable motion slightly moves towards the stable stationary nontrivial branch (see Figure 11(d)). As the fiber orientation angle increases from 4
° to 9
°, 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 (1
0
°
-
1
1
°), the probable motion goes back towards the stable stationary nontrivial branch (see Figures 11(k)-11(l)). When θ
1
=
1
2
°, 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 θ
1
=
1
6
° (see Figure 11(q)).Stationary joint probability density of the system to different fiber orientation anglesθ
1 of the first lamina with respect to the x-axis of the beam. Ω
=
5.4200, f
0
ω
2
=
3.2642, and γ
=
0.12.
(a)
θ
1
=
0
°
(b)
θ
1
=
1
°
(c)
θ
1
=
2
°
(d)
θ
1
=
3
°
(e)
θ
1
=
4
°
(f)
θ
1
=
5
°
(g)
θ
1
=
6
°
(h)
θ
1
=
7
°
(i)
θ
1
=
8
°
(j)
θ
1
=
9
°
(k)
θ
1
=
1
0
°
(l)
θ
1
=
1
1
°
(m)
θ
1
=
1
2
°
(n)
θ
1
=
1
3
°
(o)
θ
1
=
1
4
°
(p)
θ
1
=
1
5
°
(q)
θ
1
=
1
6
°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. Figures12(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 x and y in both Figures 11(f) and 11(p) to the amplitude a, 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θ
1 of the first lamina with respect to the x-axis of the beam, Ω
=
5.4200, f
0
ω
2
=
3.2642, and γ
=
0.12.
(a)
θ
1
=
5
°
(b)
θ
1
=
1
5
°


## 2.2.1. Stability Analysis of the Trivial Response

It can be found that (28a) and (28b) have a trivial steady state solution a
=
0. In order to determine the stability of the trivial steady state response, the trivial response is expanded near a
=
0 and the corresponding nonlinear terms are neglected; the linearization form of (28a) and (28b) can finally be expressed as(29a)
D
1
a
=
-
ς
a
+
f
~
4
a
sin
β
,
(29b)
D
1
β
=
σ
+
f
~
2
cos
⁡
β
+
γ
-
D
1
W
.Letr
=
ln
⁡
a and we rewrite (29a) and (29b) in the form of Itô equation as given below(30a)
d
r
=
(
-
ς
+
f
~
4
sin
β
)
d
T
1
,
(30b)
d
β
=
(
σ
+
f
~
2
cos
⁡
β
)
d
T
1
+
γ
-
d
W
.It can be found that the stochastic processβ generated on [
0,2
π
] by (30a) and (30b) is Markov and is also ergodic on [
0,2
π
] since the diffusion process is nonsingular. Naturally, the invariant measure, that is, the steady state probability density function p
(
β
) of the process β, can be determined by the following FPK equation:
(31)
d
2
p
(
β
)
d
β
2
-
d
d
β
[
(
σ
-
+
f
-
cos
⁡
β
)
p
(
β
)
]
=
0
,
where σ
-
=
2
σ
/
γ
-
2 and f
-
=
f
~
/
γ
-
2.The unique solution of (31) satisfies both the periodicity condition p
(
β
)
=
p
(
β
+
2
π
) and the normality condition ∫
0
2
π
p
(
β
)
d
β
=
1, respectively. Thus, the solution of (31) can be expressed as
(32)
p
(
β
)
=
1
C
e
(
σ
-
β
+
f
-
sin
β
)
∫
β
β
+
2
π
e
[
-
(
σ
-
x
+
f
-
sin
x
)
]
d
x
,
where C
=
4
π
2
e
(
-
σ
-
π
)
I
i
σ
-
(
-
f
-
)
I
-
i
σ
-
(
-
f
-
) and I
n
(
x
) 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
(33)
λ
(
a
0
,
β
0
)
=
lim
⁡
T
1
→
∞
1
T
1
ln
⁡
|
a
(
T
;
a
0
,
β
0
)
|
w
.
p
.
1
,
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 λ
=
λ
max
⁡; that is, λ
=
0 is the bifurcation point of the stability of the trivial response. According to [30], the largest Lyapunov exponent λ is given as
(34)
λ
=
lim
⁡
T
1
→
∞
1
T
1
ln
⁡
|
a
(
T
1
)
a
(
0
)
|
=
-
ς
+
lim
⁡
T
1
→
∞
1
T
1
∫
0
T
1
f
~
4
sin
β
(
τ
)
d
τ
=
-
ς
+
f
~
4
E
[
sin
β
]
=
-
ς
+
f
~
4
∫
0
2
π
sin
β
p
(
β
)
d
β
=
-
ς
+
f
~
4
[
I
1
-
i
σ
-
(
f
-
)
I
-
i
σ
-
(
f
-
)
+
I
1
+
i
σ
-
(
f
-
)
I
i
σ
-
(
f
-
)
]
.The variation of the largest Lyapunov exponentλ as λ
-
(
B
,
f
0
) 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.01, 0.1, and 1 and ε
ς
=
0.1, 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 B
=
2, 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 λ
=
0 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 f
0 but intensify the instability of the trivial response for a higher amplitude f
0. Such results can also be found in Figure 5 for three isohypse curves of λ
=
0 resulting from different γ; that is, γ
=
0.01, γ
=
0.1, and γ
=
1.0, respectively.The largest Lyapunov exponentλ of the trivial response of the system: γ
=
0.01, ε
ς
=
0.1. (a) Mesh surface of λ
-
(
B
,
f
0
); (b) isohypse curves of λ.
(a)
(b)The largest Lyapunov exponentλ of the trivial response of the system: γ
=
0.1, ε
ς
=
0.1. (a) Mesh surface of λ
-
(
B
,
f
0
); (b) isohypse curves of λ.
(a)
(b)The largest Lyapunov exponentλ of the trivial response of the system: γ
=
1.0, ε
ς
=
0.1. (a) Mesh surface of λ
-
(
B
,
f
0
); (b) isohypse curves of λ.
(a)
(b)Figure 5
Isohypse curves ofλ
=
0.


## 2.2.2. Stochastic Jump and Bifurcation

As the reduced (28b) has the coupled term a
D
1
β, 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]
(35)
x
=
a
sin
β
,
y
=
a
cos
⁡
β
,
and, introducing it into (28a) and (28b), we have(36a)
D
1
x
=
-
ς
x
-
f
~
4
x
y
x
2
+
y
2
-
σ
y
+
1
4
α
e
y
(
x
2
+
y
2
)
-
y
γ
-
D
1
W
,
(36b)
D
1
y
=
-
ς
y
+
f
~
4
(
y
2
+
2
x
2
)
x
2
+
y
2
+
σ
x
-
1
4
α
e
x
(
x
2
+
y
2
)
+
x
γ
-
D
1
W
,where x and y are the two perpendicular components of a; that is, a
2
=
x
2
+
y
2.Equations (36a) and (36b) are equivalent to the following set of Itô equation:(37a)
d
x
=
m
1
(
x
,
y
,
T
1
)
d
T
1
+
σ
11
(
x
,
y
,
T
1
)
d
W
+
σ
12
(
x
,
y
,
T
1
)
d
W
,
(37b)
d
y
=
m
2
(
x
,
y
,
T
1
)
d
T
1
+
σ
21
(
x
,
y
,
T
1
)
d
W
+
σ
22
(
x
,
y
,
T
1
)
d
W
,where
(38)
m
1
(
x
,
y
,
T
1
)
=
-
ς
x
-
f
~
4
x
y
x
2
+
y
2
-
σ
y
+
1
4
α
e
y
(
x
2
+
y
2
)
-
1
2
γ
-
2
x
,
b
11
(
x
,
y
,
T
1
)
=
y
2
γ
-
2
,
b
22
(
x
,
y
,
T
1
)
=
x
2
γ
-
2
,
m
2
(
x
,
y
,
T
1
)
=
-
ς
y
+
f
~
4
(
y
2
+
2
x
2
)
x
2
+
y
2
+
σ
x
-
1
4
α
e
x
(
x
2
+
y
2
)
-
1
2
γ
-
2
y
,
b
12
(
x
,
y
,
T
1
)
=
b
21
(
x
,
y
,
T
1
)
=
-
x
y
γ
-
2
,
b
i
j
(
x
,
y
,
T
1
)
=
∑
l
=
1
2
σ
i
l
(
X
,
Y
,
T
1
)
σ
j
l
(
X
,
Y
,
T
1
)
|
X
=
x
Y
=
y
,
h
=
h
h
h
h
h
h
h
h
h
h
h
h
h
h
i
,
j
=
1,2
,
and the last terms in m
1 and m
2 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
(39)
∂
p
∂
T
1
=
-
∂
∂
x
[
m
1
p
]
-
∂
∂
y
[
m
2
p
]
+
1
2
∂
2
(
b
11
p
)
∂
x
2
+
1
2
∂
2
(
b
22
p
)
∂
y
2
+
1
2
∂
2
(
b
12
p
)
∂
x
∂
y
+
1
2
∂
2
(
b
21
p
)
∂
y
∂
x
,
where p
=
p
(
x
,
y
,
T
1
∣
x
0
,
y
0
) is the transition probability density of x and y.The initial condition of FPK equation (39) is
(40)
p
=
δ
(
x
-
x
0
)
δ
(
y
-
y
0
)
,
T
1
=
0
.
The boundary conditions with respect to x and y are
(41)
p
,
∂
p
∂
x
,
∂
p
∂
y
⟶
0
,
as
|
x
|
,
|
y
|
⟶
∞
.Whenγ
-
=
0 and φ
=
0, (28a) and (28b) is a deterministically excited system and it has one trivial solution and two possible steady state nontrivial solutions when f
~
>
4
ς as given below
(42)
a
0
=
4
(
B
-
2
)
3
ε
α
±
2
(
ε
f
~
)
2
-
(
4
ε
ς
)
2
3
ε
α
.The theoretical frequency-response curves based on (42) are shown in Figure 6 for the aforementioned beam when γ
-
=
0 and φ
=
0. 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 whenf
0
=
1 and ε
ς
=
0.1.Figure7 lists a series of change of the stationary joint probability density of x and y, which is numerically solved from the FPK equation (39) by using finite difference method [34, 35] to different excitation central frequency B when γ
=
0.12 and f
0
=
1. Numerical investigation shows that the numerical result becomes convergent when the grid number is greater than 65
×
65 within the amplitude region [
(
-
2
,
2
)
,
(
-
2
,
2
)
]; thus we choose the grid number M
×
N
=
71
×
71 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 B 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 B 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 B 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 frequenciesB: γ
=
0.12, f
0
=
1. (a) B
=
2.3; (b) B
=
3.0; (c) B
=
3.5; (d) B
=
4.0; (e) B
=
5.0; (f) B
=
6.0.
(a)
B
=
2.3
(b)
B
=
3.0
(c)
B
=
3.5
(d)
B
=
4.0
(e)
B
=
5.0
(f)
B
=
6.0Figure8 reveals the stationary joint probability density to different bandwidths when B
=
2.8 and f
0
=
1.0. Hear the choice of B
=
2.8 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, γ
=
0.12 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)).Stationary joint probability density of the system to different bandwidthsγ. B
=
2.8, f
0
=
1.
(a)
γ
=
0.12
(b)
γ
=
0.20
(c)
γ
=
0.30
(d)
γ
=
0.50The 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 whenB
=
2.7.Figure10 shows the stationary joint probability density of the system to different excitation forces when B
=
2.7 and γ
=
0.10. We can find that there is only an outer flabellate peak when the excitation force, for instance, f
0
=
2.0, 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.B
=
2.7, γ
=
0.10.
(a)
f
0
=
2.0
(b)
f
0
=
1.4
(c)
f
0
=
1.0
(d)
f
0
=
0.6In order to bring forth the influence of the layup or laminate properties, in what follows, the fiber orientation angleθ
1 of the first lamina with respect to the x-axis of the beam is taken into consideration. Table 1 lists a series of data for the dimensionless natural frequency ω, the dimensionless excitation amplitude f
0, the dimensionless central frequency B, and the nonlinear parameter ε
α with a smaller change of the fiber orientation angle θ
1. 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 f
0
ω
2, and the bandwidth γ in (18) are fixed; here we choose Ω
=
5.4200, f
0
ω
2
=
3.2642, and γ
=
0.12.Table 1
System and excitation parameters corresponding toθ
1.
θ
1
ω
α
^
ε
α
f
0
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.2324Figure11 shows the stationary joint probability density of the system to different fiber orientation angles when Ω
=
5.4200, f
0
ω
2
=
3.2642, and γ
=
0.12. 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 θ
1
≤
2
° (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 θ
1
=
3
°, the probable motion slightly moves towards the stable stationary nontrivial branch (see Figure 11(d)). As the fiber orientation angle increases from 4
° to 9
°, 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 (1
0
°
-
1
1
°), the probable motion goes back towards the stable stationary nontrivial branch (see Figures 11(k)-11(l)). When θ
1
=
1
2
°, 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 θ
1
=
1
6
° (see Figure 11(q)).Stationary joint probability density of the system to different fiber orientation anglesθ
1 of the first lamina with respect to the x-axis of the beam. Ω
=
5.4200, f
0
ω
2
=
3.2642, and γ
=
0.12.
(a)
θ
1
=
0
°
(b)
θ
1
=
1
°
(c)
θ
1
=
2
°
(d)
θ
1
=
3
°
(e)
θ
1
=
4
°
(f)
θ
1
=
5
°
(g)
θ
1
=
6
°
(h)
θ
1
=
7
°
(i)
θ
1
=
8
°
(j)
θ
1
=
9
°
(k)
θ
1
=
1
0
°
(l)
θ
1
=
1
1
°
(m)
θ
1
=
1
2
°
(n)
θ
1
=
1
3
°
(o)
θ
1
=
1
4
°
(p)
θ
1
=
1
5
°
(q)
θ
1
=
1
6
°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. Figures12(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 x and y in both Figures 11(f) and 11(p) to the amplitude a, 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θ
1 of the first lamina with respect to the x-axis of the beam, Ω
=
5.4200, f
0
ω
2
=
3.2642, and γ
=
0.12.
(a)
θ
1
=
5
°
(b)
θ
1
=
1
5
°


## 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θ
1 of the first lamina with respect to the x-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 θ
1 increases from 0
° to 1
6
°.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.



---

*Source: 617828-2014-04-14.xml*