Accepted Manuscript

Advanced Nonlinear Dynamic Modelling of Bi-stable Composite Plates

Zhangming Wu, Hao Li, Michael I. Friswell

PII: S0263-8223(18)31034-1
DOI: https://doi.org/10.1016/j.compstruct.2018.06.072
Reference: COST 9866

To appear in: Composite Structures

Received Date: 16 March 2018

Revised Date: 1 June 2018
Accepted Date: 16 June 2018

Please cite this article as: Wu, Z., Li, H., Friswell, M.I., Advanced Nonlinear Dynamic Modelling of Bi-stable
Composite Plates, Composite Structures (2018), doi: https://doi.org/10.1016/j.compstruct.2018.06.072

This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers
we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and
review of the resulting proof before it is published in its final form. Please note that during the production process
errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.
 Advanced Nonlinear Dynamic Modelling of Bi-stable
Composite Plates
Zhangming Wu1,2*, Hao Li3*, Michael I. Friswell4*
1
Cardiff School of Engineering, Queens Buildings, The Parade, Newport Road,
Cardiff CF24 3AA, UK
2
School of Aerospace Engineering and Applied Mechanics, Tongji University, 1239
Siping Road, Shanghai 200092, P. R. China
3
Shanghai Institute of Satellite Engineering, P. R. China
4
College of Engineering, Swansea University Bay Campus, Fabian Way, Crymlyn
Burrows, Swansea SA1 8EN, UK
Corresponding Authors: Zhangming Wu (wuz12@cardiff.ac.uk) and
Hao Li (hithaoli@126.com) and Michael I. Friswell (m.i.friswell@swansea.ac.uk)

Abstract
This paper proposes a novel analytical model to study the nonlinear dynamics of
cross-ply bi-stable composite plates. Based on Hamilton’s principle, in conjunction
with the Rayleigh-Ritz method, an advanced analytical model with only 17 unknown
terms is developed to predict the entire nonlinear dynamic response of bi-stable
composite plates, which are excited by an electrodynamic shaker. The coupling
between the bi-stable plate and the shaker is considered in the development of
analytical model. This work, for the first time, simulates the full dynamics of bi-stable
plates using an analytical model, including the prediction of the nonlinear
characteristics of single well vibration and cross well vibration. Numerical results on
three vibrational patterns of two standard cross-ply composite plates are obtained to
study the nonlinear dynamics of bi-stable plates. The prediction accuracy on the
dynamic characteristics of different vibrational patterns of bi-stable plates are verified
by both finite element analysis (FEA) and experimental results. Large amplitude
cross-well vibrations due to the transitions between different stable states of bi-stable
plates are also characterized accurately. Applying this 17-term analytical model for
the dynamic analysis of bi-stable plates is straightforward, as the mass and stiffness
properties are obtained directly from the geometry and material properties. Only the
damping coefficients for different plates need to be determined from experiments.
Furthermore, this proposed 17-term analytical model has much higher computational
efficiency than FEA.

1. Introduction

Thin unsymmetric composite plates which exhibit multiple stable static

configurations have attracted considerable attention for over three decades [1-4]. A
bi-stable fiber reinforced composite plate possesses two stable equilibria, and the plate
is able to maintain either of the equilibria without external forces. The bi-stability of
composite plates enables large deflections from one equilibrium state to another with
small energy input. The static characteristics of unsymmetric bi-stable fiber reinforced

1
plates have been well quantified, after a comprehensive research effort over the last
three decades. Bi-stable composite plates have shown extensive potential applications
for morphing aircraft [5-10] and energy harvesting devices [11-13].
Morphing structures change shape to continually optimize operating conditions,
according to the environment or loads. In the morphing concepts that take advantage
of bi-stable plates, the equilibrium states are matched to shape requirements. Several
aero-structures applications with bi-stable morphing airfoils have been proposed and
have been experimentally and numerically investigated [5-10]. More recently,
research works shed the lights on the dynamic excitation that can trigger the
bi-stability of composite plates from one equilibrium state to another with small
energy input, for example [14-17]. Recently, the dramatic increase of wireless sensors
and electronics, has meant that energy harvesting from ambient vibration has become
an active research topic [18, 19]. Nonlinear piezoelectric harvesters have been
proposed for broadband energy harvesting. Energy harvesters that employ bi-stable
plates as the supporting structure have been proposed. The natural characteristics of
bi-stable plates, including the nonlinear bending stiffness and the snap-through
phenomenon, are used to increase the power generated and increase the effective
frequency range of operation. In particular, the snap-through between equilibrium
states gives large velocities that increases the output power [19]. For morphing
applications, bi-stable plates are inevitably exposed to high levels of dynamic
excitation, and early failure may be induced. Moreover, undesired snap-through may
be triggered by the dynamic excitation. For energy harvesting applications, complex
vibration patterns including single-well oscillations, chaos and limit cycles are
observed in experiments [20]. Lots of research interests had been attracted in the last
decade on the investigation of nonlinear dynamic characteristics of bi-stable structures,
due to the potential applications in energy harvesting [18-21]. However, modeling the
complex nonlinear dynamics of bi-stable composite plates is still an open challenge.
The purpose of this article is to derive a reduced-order analytical model to predict
the complex nonlinear dynamics of cross-ply bi-stable composite plates. The
nonlinear dynamics of bi-stable plates confined to a single stable state has been
experimentally studied by Arrieta et al. [22]. A simple model based on the observed
results of two points is established, and captures the complex dynamic phenomenon
of the tested bi-stable plates. However the model cannot predict the cross-well
dynamics. Later, Arrieta et al. [23] extended the simple model in Ref. [22] by
including a piece-wise restoring force in the nonlinear modal equations, to account for
the cross-well dynamics of bi-stable plates. It is important to understand the transient
shape of bi-stable plates during cross-well oscillations, and the simple models in Refs.
[22] and [23] cannot provide the overall response of the plate. Diaconu et al. [24]
established an analytical dynamic model of a bi-stable plate using Hamilton’s
principle, and used it to study the dynamics of the snap-through. The oscillation after
snap-through is qualitatively predicted by the model, however, it over predicts the
snap-through load by 30% compared to FEA and experiments. Based on Hamilton’s
principle, Vogl and Hyer [25] established the linear mode prediction model of
bi-stable plates, by only retaining the linear terms. The developed analysis generally
captures the major features of the lowest vibration mode of the plates. Arrieta et al.
[26] developed a low order model for the dynamics of cross-ply plates confined to a
stable state. The mode shapes of the stable state are obtained by a Galerkin approach
projecting the solution of the nonlinear problem onto the mode shapes of the linear
problem. However, the parameters corresponding to the nonlinear terms of the model

2
had to be identified from experimental frequency response functions. Taki et al. [27]
established an analytical model to analyze the nonlinear dynamic behavior of bi-stable
composite laminates with [0–90]T and piezoelectric layers using the Rayleigh–Ritz
method and Hamilton’s principle. In their modelling, a fourth-order polynomial set of
shape functions was used to achieve improved accurate results, which are validated
with the finite element method.
A desired dynamic model of a bi-stable plate should be able to quantitatively
predict the overall response of its complex dynamics, provided the dimensions,
material properties, environmental parameters and the excitation are provided. This
study establishes a dynamic model for cross-ply bi-stable plates with high
computational efficiency and high accuracy, directly from the plate geometry and
material properties. A 17 degrees of freedom dynamic model, derived using
Hamilton’s principle, is used to investigate the dynamics of bi-stable plates excited at
plate center in the transverse direction. Different vibration patterns, distinguished by
the snap through phenomenon, are predicted and the corresponding characteristics are
discussed. The accuracy of the proposed model is validated by FEA and experiment.

2. Model Derivation

Hamilton’s principle is employed to establish the dynamic model for cross-ply

bi-stable composite plates. The variation of the Lagrangian function, integrated with
respect to time, is set to zero. Thus
t2
∫t δ (T ( t ) + WF ( t ) − Π ( t ) ) dt = 0
1
(1)

where T ( t ) is the kinetic energy, WF ( t ) is the work done by the applied force, and
Π ( t ) is the total potential energy.

A typical cross-ply bi-stable plate has two stable cylindrical states and one
unstable saddle state. The x-y-z Cartesian coordinate system is employed to describe
the problem, as illustrated in Fig. 1. The origin of the coordinate system is at the
center of the plate. The x and y axes are parallel or perpendicular to the fiber
directions.

Fig. 1. Sketch of the unstable and stable configurations of a bi-stable cross-ply composite plate
and the analysis coordinate system (reproduced from Hyer [3])

The bi-stable plates are assumed to be thin. Thus, the Kirchhoff hypothesis is
employed and the layers are assumed to be in the state of plane stress [28]. Hence

3
  ε  ε x0   k 
 x   0   x 
 ε y  = ε y  + z  k y  (2)
   0  
γ xy  γ xy  k xy 

where

 ∂2w 
− 2 
 k   ∂x 
 x   ∂ 2 w 
k =  ky  =  − 2  (3)
   ∂y 
 k xy   2 
−2 ∂ w 
 ∂x∂y 
represents the plate curvatures. The deflections of bi-stable plates are relative large
compared to the plate thicknesses, and hence the geometrical nonlinearity has to be
considered in the geometric equation [29]. The mid-plane strains, ε 0 , are then
defined in the spirit of von Karman, and the in-plane strains of a thin bi-stable plate
take the form
 ∂u 0 1  ∂w 2 
 +   
 0   ∂x 2  ∂x  
 εx   2 
0  0   ∂v
0
1  ∂w  
ε = εy  =  +    (4)
 0   ∂y 2  ∂y  
γ   0 0 
 xy   ∂u + ∂v + ∂w ∂w 
 ∂y ∂x ∂x ∂y 
 

where u0, v0 and w denote the in-plane and out-of-plane displacements. The total
potential energy of the plate under a thermal-induced load is given by
Lx 2 Ly 2  1  0T  A B  ε 0   T ε 0  
kT  
Π=
∫ ∫
− Lx 2 − Ly 2
2 ε
 

 − N
  B D  k   s
M Ts     dydx
  k 
 
(5)

where Lx and L y denote plate length and width, respectively, and the explicit time
dependence has been omitted. A, B and D are the stretching stiffness matrix, the
stretching-bending coupling matrix and the bending stiffness matrix, respectively.
N s and Ms are the resultant force and moment due to the thermal stress. The
0
potential energy depends on the in-plane strains ε , and the curvatures k, which are
determined by the assumed displacement functions, u0, v0 and w.
0
The mid-plane strain ε , and the bending curvature, k, for an unsymmetric plate
are coupled, since the stretching-bending coupling matrix, B, is non-zero. If a

4
displacement model for the out-of-plane displacement function, w, is assumed, then
the mid-plane strains, ε x0 and ε 0y , should have terms which are dependent on w, due
to the stretching-bending coupling. However, previous studies often applied a
generalized modeling approach for unsymmetric plates, where ε x0 , ε 0y and w are
defined by independent series [30, 31]. However, the influence of stretching-bending
coupling on the assumed form of the mid-plane strains, ε x0 and ε 0y , were not
explicitly considered. As a result, more terms have to be used in the series expansion
of mid-plane strains to achieve better prediction accuracy.
In this study, an alternative procedure is applied to model the mid-plane strains
0
ε for cross-ply bi-stable plates. According to classical lamination theory, the
in-plane stress resultants for a thermally-loaded unsymmetric plate is expressed as

N = A ε 0 + B k − N Ts (6)

Equation (6) can be rearranged to give

ε 0 = A -1N + A -1N Ts − A-1B k = ε m + A-1N Ts − A -1B k (7)

where A-1N can be essentially considered as a membrane strain field εm .

Substituting Eq. (7) into Eq. (5), the total strain energy becomes
Lx 2 Ly 2  1  mT A 0  ε m 
kT  
Π=
∫ ∫
− Lx 2 − Ly 2
2 ε
 

 0 D − BA -1B   k 
  
(8)
 1 A -1 − A -1B   N s  
−  N Ts M Ts   2     dydx
    k  
 0 I

This process of reformulating the strain energy into Eq. (8) is analogous to the
reduced stiffness matrix (RSB) method [32]. In doing so, the elastic component of
1 mT
total strain energy is decoupled into a stretching part ε Aε m and a bending part
2
1 
k D − BA -1B  k . Therefore, the deformation of a plate can be decoupled into
2  
independent stretching and bending, while the neutral axis does not necessarily lie at
the mid-plane, as for a standard isotropic thin shell.
Models that assume constant bending curvatures over predict the plate stiffness
[24], and hence non-constant curvatures are assumed in this study. The sixth order
out-of-plane displacement w( t ) is defined, and the unknown coefficients are again
time dependent for the dynamic problem. Thus

w ( t ) = a ( t ) x 2 + b ( t ) y 2 + a1 ( t ) x 4 + b1 ( t ) y 4 + a2 ( t ) x6 + b2 ( t ) y 6 + e ( t ) x 2 y 2 (9)

where the unknown coefficients, ai ( t ) , bi ( t ) and e ( t ) , represent generalized

co-ordinates in the model. The form of w( t ) is similar with the form proposed by
Gigliotti et al. [33], although the extra term e ( t ) x 2 y 2 is included to account for the

5
anticlastic phenomenon near the plate corners.

In this work, instead of taking the mid-plane strains ε0 as the primary variables
in the model, the dynamic stretching strains εm are assumed and expanded up to a
sixth order polynomials as

ε xm ( t ) = c ( t ) + c1 ( t ) y 2 + c2 ( t ) y 4 + c3 ( t ) y6 + c4 ( t ) x2 y 2
2 4 6 2 2
(10)
εm
y ( t ) = d ( t ) + d1 ( t ) x + d2 ( t ) x + d3 ( t ) x + d4 ( t ) x y

where the unknown coefficients, ci ( t ) and di ( t ) , represent generalized

co-ordinates in the model. Consequently, the nonlinear dynamic problem of the
cross-ply composite plate is modeled by the three variables given in Eqs. (9) and (10),
in which there are 17 unknown terms in total.
−1 T
For [0n/90n]T cross-ply composite plates, A −1B and A N s are explicitly
expressed as,

 B*11 B*12 0   B*11 B*12 0

   
− A -1B = B −1 =  B*21 B*22 0  =  − B*12 − B*11 0 (11)
   
 0 0 0  0
 
0 0


 A−1 N sT1 
−1 T  −1 T −1 T 
A Ν s =  A N s 2 = A N s1  (12)
 A−1 N sT3 = 0 
 
Substituting Eqs. (11) and (12) into Eq. (7), the mid plane strains are given by

ε x0 = ε xm + B*11k x +B*12k y + A-1N sT1

ε 0y = ε ym − B*12 kx − B*11k y + A-1 NsT1 (13)
0 m
γ xy = γ xy

0 0
The mid plane displacement functions, u ( t ) and v ( t ) , are then derived by
integrating with x and y, as

 1  ∂w ( t )  
2

∫
u (t ) =  ε x (t ) −    dx
0 0
 2  ∂x  
 
(14)
 1  ∂w ( t )  
2

∫
0 0
v (t ) = ε y (t ) −   dy
 2  ∂y  
 
m
The shear component of mid plane and stretching strains γ xy ( t ) , γ xy0 (t ) , are derived
0 0
by substituting the series forms of u ( t ) , v ( t ) and w( t ) into Eq. (2). The

6
 0 0
polynomial expressions for u ( t ) , v ( t ) and γ xy
0
( t ) are given in the Appendix.
0 0
Note, the constants in u ( t ) and v ( t ) due to integration of strains are treated as 0.

Substituting the expressions for ε xm ( t ) , ε ym ( t ) and γ xy

m
( t ) into Eq. (8), the
dynamic potential energy Π ( t ) is derived. Since the mid plane shear strain γ xy
0
(t )
is derived explicitly from the displacement functions, the compatibility of condition
for the mid-plane strains is satisfied automatically.
The bi-stable composite plates in this study are mounted and excited at the central
point. The kinetic energy is defined as
Lx Ly 2
1  d (w0 (t ) + w(t )) 
T (t ) = ρ h ∫ Lx ∫ 2Ly 
2
 dydx (15)
2 2  dt 
− −
2

where ρ denotes the mass density, w0 ( t ) denotes the excitation displacement at the
central point, and h is the thickness of the plate. In Eq. (15), the energy contributions
given by the transverse inertia terms are neglected.
In the present study, four constant and equal concentrated transverse forces are
applied at the four corners of the plate. In this case, the work done by the applied
forces, WF, is
WF ( t ) =
 2
 Ly 
2 4
 Ly 
4 6
 Ly 
6 2 2
  Lx   Lx   Lx   Lx   Ly   (16)
4F a(t)   +b(t)   +a1(t)  +b1(t)  +a2(t)  +b2(t)  +e(t)   
 2 2 2 2 2 2 2 2 
 

3. Excitation Controlled by Displacement

In the present study, a sinusoidal excitation is applied at the center of the bi-stable
plate in the transverse direction. If the bi-stable plate is connected to a structure with
great mass and stiffness, the influence from the vibration of the supporting structure
can be ignored. In this case, a prescribed excitation displacement w0 ( t ) is
introduced to represent the motion of the supporting structure in the dynamic model.
Since there are 17 degrees of freedom in the proposed dynamic model, the generalized
displacement is expressed as
T
X1 ( t ) = [ a , b , a1 , b1 , a 2 , b2 , e , c , d , c1 , d1 , c2 , d 2 , c3 , d 3 , c4 , d 4 ] (17)

According to the Hamilton principle, the variation of the Lagrangian function,

given by (T − Π − WF ) , is equal to zero. The 17 equations of motion are thus derived,
and expressed in the following matrix form

7


+ D ( X
M1X
)+ K (X ) = F (18)
1 1 1 1 1

where M1 is the mass matrix, D ( X

) is the damping force, K ( X ) is the
1 1 1
nonlinear stiffness force, and the overdot denotes differentiation with respect to time.
F1 denotes the excitation force and also can be expressed by F1 = Gw

0 , which is the
base excitation force, which depends on the plate inertia and the acceleration w

. The
vector G is given in the appendix.
With the definition of the generalized co-ordinates in Eq. (17), and the definition of
0 0
the displacement models for u ( t ) , v ( t ) and w ( t ) given by Eqs. (9) and (14),
the mass matrix M1 is easily derived. Similarly, the stiffness matrix K 1 is derived by
writing the strain energy of the plate, Π , in terms of X1 , and differentiating in the
usual way. The expressions of the matrices M1 and K 1 are given in Appendix.
Since the displacement w ( t ) is a linear function of the degrees of freedom, the mass
matrix M1 is a constant coefficients matrix.

Damping is very difficult to model precisely, and here the assumption of Rayleigh
proportional damping [24] is employed and thus:

D ( X1 , X

) =α M X
1


1 1 (19)

where α is the mass proportional damping coefficient. As the vibration of bi-stable

plates often occurs at low frequencies, the contribution of stiffness damping is ignored.
Thus, only the mass damping effect is considered in this study, and the damping force
is a function of velocity alone.

4. Excitation Controlled by Force from a Shaker

In previous studies on the dynamics of bi-stable plates, electrodynamic shakers

were often used to excite the plates. Arrieta et al. used a transducer between the plate
and the shaker to measure the excitation force, and a linear feedback proportional
controller was implemented to eliminate the dynamic coupling between the plate and
the shaker [23, 34]. As the dynamics of bi-stable plates is nonlinear and complex, it is
unlikely that the coupling between the plate and the shaker can be eliminated,
completely. In this study, the coupling between the plate and the shaker is explicitly
considered using the following analytical model. An electromagnetic shaker consists
of a moving part, the armature assembly, a linear spring support for the armature and
a dashpot, as illustrated in Fig. 2.

8
 Fig. 2 Dynamic model of the electromagnetic shaker

The single degree of freedom dynamic model of the shaker is given by

0 + Cw
0 + kw0 = f ( t )
mw (20)

where m is the mass of the moving armature, C is the damping coefficient of the
dashpot, k is the stiffness of the spring system, and f ( t ) denotes the sinusoidal
electromagnetic force of the shaker. Due to the coupling between the shaker and the
plate, the oscillation amplitude of the center of the plate, w0 , is unknown. Therefore,
the system is a combination of the bi-stable plate and the shaker. There are now 18
degrees of freedom as w0 is now also included. Thus,

 w0 
X2 =   (21)
 X1 
The dynamic equations are derived similarly as the process of deriving Eq. (18),
M2X
+ K (X ) = F


+ CX (22)
2 2 2 2 2

Each term in Eq. (22) is expressed in matrix form as

C 0 
C=   (23)
 0 α M1 

 kX 2 (1) 
K 2 ( X2 ) =   (24)
K1 ( X1 ) 
T
F2 ( t ) =  f ( t ) 0 ... 0 (25)

the complete form of M2 is presented in the Appendix

9
5. The Effect of Inertial Mass

Inertial masses are often attached to the corners of bi-stable plates to increase the
dynamic response. The contribution of these masses (assumed to be equal at each
corner) to the kinetic energy is

1  4 2 4 2 4 2
Tinertia = minertial  ∑ w
i + ∑ u
i + ∑ v
i  (26)
2  i =1 i =1 i =1 

where the velocities of each mass (i.e. each corner) are

w
1,…,4 ( t ) = w
( ±
1
2
1
)
Lx , ± Ly , t + w
0 ( t )
2

u
1,…,4 ( t ) = u
(
0
±
1
2
Lx , ±
1
2
L ,t )
y (27)

v
1,…,4 ( t ) = v
(
0
±
1
2
Lx , ±
1
2
L ,t)
y

and minertial represents the mass added to each corner of the plate. The contribution to
the system mass matrix is easily derived. In Eq. (18), the contribution of inertial mass
is added to the mass matrix by replacing M1 with M1+Minertial1. Similarly, the mass
matrix M2 should be replaced by M2+Minertial2 in Eq. (22). The matrices Minertial1 and
Minertial2 are given in the Appendix. For the excitation force in Eq. (18), the
contribution of the inertial mass F ( X1 )inertial is given in the Appendix, and F ( X1 )
in Eq.(8) is replaced by F ( X1 ) + F ( X1 )inertial .

6. Experimental Setup

In the present study, an experiment is established to test the dynamic response of

bi-stable composite plates. Figure 3(a) illustrates the schematic diagram of this
experimental system, in which each experimental device is connected with respect to
the signal (data) flow. Fig. 3(b) shows the practical experimental set up, in which the
bi-stable plate is mounted at its centre and excited by a shaker. An accelerometer is
attached at the centre of the bi-stable plate to measure the excitation acceleration. A
laser vibrometer is used to measure the full dynamic response of the bi-stable plate,
and particularly the motion at the plate corners where the 11.54 g inertial masses are
attached.

10
 Data
collecter

Inertial mass

Bi-stable plate

PC Acceleration
transducer

shaker

Laser
vibrometer Shaker
controller

(a) Schematic diagram for the experimental set up

(b) A bi-stable plate with four inertial masses and mounted at the centre.
Fig. 3 The experimental setup

7. Finite Element Analysis

In this study, finite element analysis (FEA) for the bi-stable plates was performed
using the commercial tool ABAQUS to validate the proposed analytical model. The
S4R shell element (4-node general-purpose shell, reduced integration with hourglass
control, finite membrane strains) was chosen to model the bi-stable plates. A mesh
density of 40 × 40 was chosen in the finite element model (FEM) to achieve the
required accuracy and efficiency.
A static analysis with two different conditions predicted the stable configurations
for the bi-stable plates. In the first “Static” step, a stable configuration is determined

11
under the conditions that the plate cools down from the curing temperature to room
temperature, and a concentrated force along transverse direction is applied at each
corner. In the second “Static” step, another stable configuration is obtained for the
bi-stable plate, where the transverse concentrated forces are removed, the temperature
is kept constant and the second stable configuration is obtained.
After the static analysis, the nonlinear dynamic analysis is performed for the
bi-stable plates using the “Dynamic implicit” step in ABAQUS. The plate is mounted
at the centre and excited at the centre in the transverse direction. For the case of
prescribed excitation displacement, a sinusoidal displacement is applied to the centre
of the plate. For the case of excitation by the shaker, a spring and dashpot is
connected to the centre of the plate, and a sinusoidal force is applied to the plate
centre in the transverse direction, as illustrated in Fig. 4.

Fig. 4 Schematic of the finite element model of the plate and shaker system

8. Results and Discussion

The bi-stable plates are manufactured using the unidirectional carbon fiber-epoxy
matrix prepreg, CCF300/5428. The material properties of the CCF300/5428 prepreg
are given in Table 1.
Table 1 Material properties of CCF300/5428 prepreg

CFRP E11=145 GPa, E22=9.75 GPa, G12=5.69 GPa, ν12 =0.312,

α11= 0.4×10-6/°C, α22 = 25×10-6/°C, Thickness=0.125mm

8.1 Static analysis

The static equilibrium shapes of the bi-stable plates need to be determined as the
initial states for dynamic analysis. The static equilibrium configurations of the
bi-stable plates will be obtained by solving the nonlinear coupled equations
K1 ( X1 ) = 0 (28)

12
The “FindRoot” function with initial values of X 1 is employed in Mathematica®
to solve Eq. (28). For the bi-stable plates, there are three groups of solutions, and each
group of solutions represents an equilibrium configuration. The actual equilibrium
configuration depends on the initial value of the unknowns X 1 . The stability of the
equilibrium solutions need to be checked by the Jacobian matrix, which is positive
definite for stable configurations:
∂ ( K1 ( X1 ) )
J= (29)
∂ ( X1 )

The static stable configurations of the square composite plates 100mm×100mm,

[02/902]T and 100mm×100mm, [03/903]T are analyzed in this section. Each plate is
mounted at its centre, and concentrated forces are applied at each corner of the
bi-stable plates. The plates cool down from the curing temperature of 180℃ to the
room temperature of 20℃. The stable configurations for the bi-stable plates are
predicted by both the analytical model and the FEA. The two stable cylindrical mode
shapes of each square bi-stable plate are identical in principle, and the only difference
lies in the different orientation of the cylindrical deformation which are perpendicular
to each other. Fig. 5 and Fig. 6 present the predicted stable states of the two bi-stable
plates [02/902]T and [03/903]T, respectively, and show that the stable configurations
predicted by the two methods coincide quite well. The numerical results for the
out-of-plane deflection at the corners of the bi-stable plates given by the analytical
model and the FEA are compared in Table 2. The small errors between the predicted
results from the two methods, as shown in Table 2, further validate the proposed
analytical model proposed in this work.

Fig. 5 Predicted stable configurations of the 100mm×100mm, [02/902]T plate. FEA results are
represented by the black dots, and the theoretical results are represented by the continuous
surfaces.

Fig. 6 Predicted stable configurations of the 100mm×100mm, [03/903]T plate. FEA results are
represented by the black dots, and the theoretical results are represented by the continuous
surfaces.

13
 Table 2 Predicted out-of-plane deflection at the corners of the bi-stable plates

Applied Force Plate Analytical (mm) FEM (mm) Error (%)

[02/902]T 6.68548 6.77379 1.30

0N
[03/903]T 3.96181 3.95916 -0.07

[02/902]T 8.12645 8.27471 1.79

0.5N
[03/903]T 4.49868 4.53591 0.82

[02/902]T 12.3557 12.4621 0.85

2N
[03/903]T 5.96985 6.09222 2.01

8.2 Free Vibration

To study the free vibration of the bi-stable plates, four identical concentrated
forces are initially applied at each corner of the plates, which are mounted at their
centres. The applied forces are the suddenly removed. For the illustration purposes,
the numerical results are only determined and presented for the 100mm × 100mm,
[03/903]T bi-stable plate. Four inertial masses of 11.54 g are attached to each corner of
the bi-stable plates to increase the amplitude of the dynamic response. The static
stable configurations of the [03/903]T bi-stable plate predicted by the analytical model
are shown in Fig. 6, and are used as the initial deformation for the dynamic analysis.
The dynamic response of the bi-stable plates requires the solution of Eq. (18). The
excitation amplitude at the centre of the plate is zero in the free vibration analysis, i.e.
‫ݓ‬଴ (‫ = )ݐ‬0. In the FEA, the concentrated forces are applied in the “Static” step, and
are removed in the “Dynamic implicit” step. Damping is not considered in either
analysis, i.e. α = 0 .
The free vibration of the [03/903]T bi-stable plate subjected to initial concentrated
forces of 2N, 0.5N and 0N are predicted. Fig. 7(a) and 7(b) present the oscillation
curves at the plate corners determined by the analytical model and the FEA,
respectively. Figure 7 shows that the analytical results and the FEA results are very
close to each other. When the initial deformation (applied force=0.5N) is small, the
plate oscillation is nearly sinusoidal. However, with larger initial deformation (applied
force=2N), the oscillation for the bi-stable plate exhibits slightly non-sinusoidal
characteristics.
The Fast Fourier Transform (FFT) of the predicted vibration presented in Fig. 7 is
shown in Fig. 8. For small amplitude excitation (0.5N) there is only one peak, which
indicates that the vibration is purely sinusoidal. The first dominant frequency for the
analytical model is 44Hz, which is slightly higher than that given by the FEA (42Hz).
With larger amplitude excitation (2N), the FFTs have two peaks. The second
dominant frequency is twice that of the first dominant frequency, and this indicates

14
the second harmonic vibration. The peak frequencies predicted by the analytical
model (40Hz and 80Hz) are slightly higher than those predicted by FEA (38Hz and
76Hz). Moreover, the FFT results in Fig. 8 also illustrate that the first dominant
frequency decreases by 4Hz when the oscillation amplitude increases (the applied
load increases from 0.5N to 2N). This indicates that the stiffness of the bi-stable plate
decreases with respect to increasing oscillation amplitude (i.e. softening), due to the
natural nonlinearity of the bi-stable plate. Although the analytical model over predicts
the dominant frequencies by about 2Hz, it accurately predicts the nonlinearity of the
bi-stable plate. If damping is considered in both the analytical model and the FEM
with α = 100 , the predicted free vibration of the 100mm×100mm, [03/903]T plate
shown in Fig. 9 is obtained. The decay rates of the vibrations predicted by the two
models are nearly identical.

0.008 0.008
Analytical-Initial force=2N FEA-Initial force=2N
0.007 Analytical-Initial force=0.5N 0.007 FEA-Initial force=0.5N
Analytical-Initial force=0N FEA-Initial force=0N
0.006 0.006
Corner deflection (m)

Corner deflection (m)

0.005 0.005

0.004 0.004

0.003 0.003

0.002 0.002

0.001 0.001

0.000 0.000
0.0 0.1 0.2 0.3 0.4 0.5 0.0 0.1 0.2 0.3 0.4 0.5
t (s) t (s)

(a) Analytical Analysis (b) FEA

Fig. 7 Predicted free vibration of the 100mm×100mm, [03/903]T plate, with zero
damping

15
 0.0005 0.0025
Analytical-Initial force=0.5N Analytical-Initial force=2N
40Hz
44Hz
0.0004 0.0020

0.0003 0.0015

Amplitude

Amplitude
0.0002 0.0010

0.0001 0.0005

80Hz
0.0000 0.0000
25 50 75 100 125 150 175 200 25 50 75 100 125 150 175 200
Frequency Frequency

(a) Analytical, initial force=0.5N (b) Analytical, initial force=2N

0.0007 0.0030
FEA-Initial foece=0.5N FEA-Initial force=2N
0.0006 42Hz 0.0025
38Hz
0.0005
0.0020
0.0004
Amplitude

Amplitude
0.0015
0.0003
0.0010
0.0002

0.0001 0.0005
76Hz
0.0000 0.0000
25 50 75 100 125 150 175 200 25 50 75 100 125 150 175 200
Frequency Frequency

(c) FEA, initial force=0.5N (d) FEA, initial force=2N

Fig. 8 FFT analysis of the predicted free vibration of the [03/903]T plate in Fig. 7

4.6 FEA
Displacement of corner node (mm)

Analytical
4.4

4.2

4.0

3.8

3.6

3.4

3.2
0.0 0.2 0.4 0.6 0.8 1.0 1.2 1.4 1.6 1.8 2.0
Time (s)

Fig. 9 Predicted free vibration of the 100mm×100mm, [03/903]T plate with the
proportional damping coefficient α = 100

8.3 Single well vibration

The forced vibration of the bi-stable plates within a single energy well is
investigated by experiment, analytical model and FEA. The experimental setup for the
vibration testing of bi-stable plates is illustrated in Fig. 3. The excitation acceleration

16
at the plate centre was monitored throughout the experiments. For the low amplitude
single well vibration, the excitation acceleration at the plate centre is sinusoidal. The
excitation amplitude is adjusted manually from zero to the target level at each
frequency. The oscillation amplitudes at the plate corner at different frequencies were
recorded by a laser vibrometer, with a fixed level of excitation acceleration for each
measurement. The measured vibration amplitudes at the plate corners for the two
bi-stable plates ([02/902]T and [03/903]T) are presented in Fig. 10.
Eq. (18) is solved to predict the small amplitude single well vibration of bi-stable
plates. Diaconu et al. [24] suggested a measured mass damping coefficient of α = 83 ,
although α depends on the plate dimensions, material properties, plate lay-up, etc.
Preliminary studies indicate that changing α at these excitation magnitudes has a
negligible effect on the stable dynamics of single well vibration if the excitation time
is sufficiently long. This benchmark study focuses on the validation of the analytical
model against the FEA rather than precisely matching the experimental results.
Therefore, a damping coefficient of α = 100 is employed. The vibration amplitudes
at the plate corners corresponding to different frequencies are predicted using the
proposed analytical model and are also shown in Fig. 10.
The single well vibrations at different frequencies are also simulated by FEA. A
sinusoidal displacement is applied at the centre of the bi-stable plates in the transverse
direction. Compared with the analytical model, the time domain analysis in FEA using
the “dynamic implicit” step is very time consuming. The oscillation amplitudes at the
plate corners for different excitation frequencies are numerically determined by FEA
and are compared with the results predicted by the analytical model, as shown in Fig.
10.
Fig. 10 shows that both the analytical model and the FEA predict only one
response peak within the frequency range, whereas extra peaks were observed in the
experimental results for each bi-stable plate. To study the origin of the extra resonant
peaks in the experimental results, a modal analysis of the bi-stable plates mounted at
the centre was performed using FEA. The mode shapes of a square [03/903]T bi-stable
plate predicted by FEA are presented in Fig. 11. For the first and second modes, the
bi-stable plate rotates around its centre. For the third mode, the plate exhibits a
symmetric bending deformation, while the plate exhibits an unsymmetric bending
deformation along the diagonal line in the fourth mode. Therefore, the common peak
in the frequency response function of the bi-stable plate that is predicted by analytical
model, FEA and experiment corresponds to the third resonance mode. In FEA, the
model is ideal and symmetric, and thus the other modes are not excited in the time
domain analysis. Since symmetric displacement function is employed in the analytical
analysis, the bi-stable plate is unable to deform in mode 1, mode 2 and mode 4, that
are shown in Fig. 11. However, in the experiment the specimen is not ideal and has
variabilities in its dimensions and symmetry, including the position of the inertial
masses. Therefore, resonance responses occur with respect to the extra vibrational
modes, other than mode 3, which were also excited and observed in the experiments.
It should be noted that the predicted amplitude-frequency trend of the single well
vibration of the bi-stable plate also exhibits nonlinear characteristics, given by a
softening stiffness characteristic with increasing vibration amplitude.

17
 0.5g,Experimental mode 3
6

Oscillation Amplitude of corner node (mm)

0.5g.Theoretical
0.5g,FEA
5
mode 4
4
mode 1&mode2
3

10 15 20 25 30 35
Frequency (Hz)

(a) 100mm×100mm, [02/902]T

Oscillation Amplitude of corner node (mm)

1.8 mode3
0.2g-Experimental
1.6 0.2g-Theoretical mode4
0.2g-FEA
1.4

1.2

1.0

0.8 mode1&mode2

0.6

0.4

0.2

0.0

30 35 40 45 50 55
Frequency (Hz)

(b) 100mm×100mm, [03/903]T

Fig. 10 Vibration amplitudes of the bi-stable plates under low level sinusoidal
excitation at different frequencies

18
 (a) initial (b) mode1 (c) mode2 (d) mode 3 (e) mode 4
shape 15.29HZ 15.61Hz 41.95Hz 46.51Hz
Fig. 11 The first four linear vibration modes about the stable configuration of the
square [03/903]T bi-stable plate

The predicted and measured displacement of the plate corner for the 100mm×
100mm, [02/902]T and 100mm×100mm, [03/903]T plates are presented in Fig. 12. The
dynamic responses predicted by the analytical model agree well with experimental
results, and the oscillation curves at the plate corner are almost perfectly sinusoidal.

Excitation frequency is 25.5Hz, excitation acceleration amplitude is 0.5g

4
Experimental
3 Analytical
FEA
corner displacement (mm)

-1

-2

-3

2.00 2.02 2.04 2.06 2.08 2.10 2.12 2.14 2.16 2.18 2.20
t (s)

(a) 100mm×100mm, [02/902]T

19
 Excitation frequency is 43Hz, excitation acceleration amplitude is 0.2g

1.0
Experimental
0.8 Theoretical
FEA

corner displacement (mm)

0.6

0.4

0.2

0.0

-0.2

-0.4

-0.6

-0.8

-1.0
2.00 2.02 2.04 2.06 2.08 2.10 2.12 2.14 2.16 2.18 2.20

t (s)

(b) 100mm×100mm, [03/903]T

Fig. 12 The corner displacement of the bi-stable plates corresponding to single-well oscillation
(dT=160℃)

8.4 Cross-well vibration

When the vibration amplitude of the bi-stable plate is larger than a critical level,
the dynamic snap-through phenomenon occurs. In other words, the bi-stable plate
transforms between the two stable configurations during the vibration. In the present
study, if the snap-through phenomenon occurs in each vibration cycle, the vibration
pattern is termed continuous cross-well vibration; otherwise, the vibration is termed
intermittent cross-well oscillation.
In the experiment, intermittent cross-well oscillation of the 100mm×100mm,
[03/903]T plate was observed. When the excitation force was increased slowly and
achieved a critical value, the vibration pattern of 100mm×100mm, [03/903]T plate
suddenly transforms from the single-well oscillation to the intermittent cross-well
oscillation. This phenomenon is clearly shown in Fig. 13, where the measured
excitation acceleration is given as the amplitude increases slowly. During the
intermittent cross-well vibration the excitation is no longer sinusoidal and has
irregular amplitude. This phenomenon indicates that the bi-stable plate and the shaker
are strongly coupled during intermittent cross-well vibration. In this section, Eq. (22)
is solved to analytically predict the cross-well vibration of bi-stable plates. In Eq. (22),
the excitation displacement w0 (t) is assumed to be unknown and the coupling
between the plate and the shaker is considered.

20
 Fig. 13 Experimental acceleration at the centre of the 100mm×100mm, [03/903]T plate, at
excitation frequency 44 Hz

The nonlinear dynamic response of the 100mm × 100mm, [03/903]T plate

predicted by the analytical model is shown in Fig. 14. The amplitude of the sinusoidal
excitation force of the shaker varies linearly with time, that is F(t) = t ×sin(44× 2πt) .
With the increase in the excitation amplitude, the vibration pattern of the bi-stable
plate transforms from a single-well oscillation to an intermittent cross-well oscillation,
and then to a continuous cross-well oscillation. The analytical model successfully
predicts the irregular pattern of the intermittent cross-well vibration at the plate centre,
which is also roughly consistent with the nonlinear characteristics observed in the
experiment. When the continuous cross-well vibration occurs, the dynamic response
of the bi-stable plate becomes regular again.

21
 F (t ) = t × Sin[44 × 2π t]
25
intermittent cross-well

Excitation acceleration of shaker (g)

20
single-well

15

10

-5

-10

-15

-20

-25
2.0 2.5 3.0 3.5 4.0 4.5 5.0 5.5 6.0 6.5 7.0 7.5 8.0

t (s)

F (t ) = t × Sin[44 × 2π t]
70
intermittent cross-well continuous
Excitation acceleration of shaker (g)

60
50 cross-well
40
30
20
10
0
-10
-20
-30
-40
-50
-60
-70
8.0 8.5 9.0 9.5 10.0 10.5 11.0 11.5 12.0

t (s)

Fig. 14 Acceleration at the centre of the 100mm×100mm, [03/903]T plate as the actuation force
amplitude linearly increases, at actuation frequency 44 Hz.（dT=160℃, C=5）

The measured and predicted corner displacements of the 100mm × 100mm,

[02/902]T and 100mm × 100mm, [03/903]T plates corresponding to intermittent
cross-well vibration are compared in Fig. 15. Although the predicted displacements do
not match exactly the experimental results, the chaotic characteristics of the
intermittent cross-well vibration are very close.

22
 10
25Hz, Analytical
5

0
Displacement of corner node (mm)
-5

-10
0.0 0.2 0.4 0.6 0.8 1.0 1.2 1.4 1.6
Time (s)

10
25Hz, Experimental
5

-5

-10
0.0 0.2 0.4 0.6 0.8 1.0 1.2 1.4 1.6
Time (s)

(a) 100mm×100mm, [02/902]T

4 44Hz, Analytical

2
Displacement of corner node (mm)

-2

-4
1.00 1.25 1.50 1.75 2.00 2.25 2.50
Time (s)

8
44Hz, Experimental
6
4
2
0
-2
-4
1.00 1.25 1.50 1.75 2.00 2.25 2.50
Time (s)

(b) 100mm×100mm, [03/903]T

Fig. 15 Displacements of the cross-ply bi-stable plates corresponding to intermittent cross-well

oscillation

23
Fig. 16 presents the predicted corner displacements of the 100mm×100mm, [02/902]T
and 100mm×100mm, [03/903]T plates during the continuous cross-well vibration.
The excitation force of the shaker is 6N. The continuous cross-well vibrations of these
two bi-stable plates also exhibit slight nonlinearity; nevertheless, the cross-well
vibration represented by the corner displacements is regular. The FFT analysis for the
cross-well vibration of these two bi-stable plates are performed and shown in Fig. 17.
The FFT results indicate the presence of the third harmonic. In the experiment, it was
difficult to observe the continuous cross-well vibration, as the deflections or
dimension asymmetry of the bi-stable plates had an adverse effect on the continuous
snap-through event.
Excitation frequency f=26Hz, excitation force amplitude=6N Excitation frequency f=40Hz, excitation force amplitude=6N
10 10
analytical analytical
Displacement of corner node (mm)

Displacement of corner node (mm)

8 FEA 8 FEA
6 6
4 4
2 2
0 0
-2 -2
-4 -4
-6 -6
-8 -8
-10 -10
0.60 0.65 0.70 0.75 0.80 0.85 0.90 0.95 1.00 0.80 0.82 0.84 0.86 0.88 0.90 0.92 0.94 0.96 0.98 1.00
t (s) t (s)

(a) 100mm×100mm, [02/902]T (b) 100mm×100mm, [03/903]T

Fig. 16 Displacement curves of cross-ply bi-stable plate corresponding to cross-well oscillation

(dT=160℃, C=5)

24
 5 2.5
26Hz
4 2.0 25Hz

3 1.5
Amplitude

Amplitude
2 1.0

1 0.5
78Hz 77Hz
0 0.0
20 40 60 80 100 120 140 160 180 200 20 40 60 80 100 120 140 160 180 200
Frequency Frequency

(a) 100mm×100mm, [02 /902]T, analytical (b) 100mm×100mm, [02/902]T, FEM

6 3.0
40Hz 40Hz
5 2.5

4 2.0
Amplitude

Amplitude
3 1.5

2 1.0

1 0.5
120Hz 120Hz
0 0.0
0 20 40 60 80 100 120 140 160 180 200 20 40 60 80 100 120 140 160 180 200
Frequency Frequency

(c) 100mm×100mm, [03 /903]T, analytical (d) 100mm×100mm, [03/903]T, FEM

Fig. 17 FFT of the corner displacements in Fig. 16

The dynamic response of the coupled system of the plate and the shaker does not
directly reflect the natural dynamics of the bi-stable plates, whereas the natural
dynamics can be obtained by exciting the uncoupled plate with a given displacement
input. Specifically, Eq. (18) is applied to study the dynamics of bi-stable plates
analytically, since the excitation acceleration of the plate is not coupled with the
shaker in the dynamic model represented by Eq. (18). The time domain responses of
the [02/902]T and [03/903]T plates excited by the sinusoidal input w0 (t ) are predicted
by Eq. (18), and the resulting vibration displacements at the plate corners are plotted
in Fig. 18. The excitation amplitude increases by 0.1g in every vibration cycle. The
excitation amplitude is linearly varying with time, and the formula for the excitation
input is also given in the caption of Fig. 18. The vibration curves in Fig. 18 clearly
show that when the excitation amplitude increases to a critical value, the vibration
pattern of both plates transfer from the linear sinusoidal vibration to the cross-well
vibration. The critical excitation acceleration amplitude, with which the vibration
pattern transforms to the cross-well oscillation, for [02/902]T plate at 22 Hz is 3.6g,
and for [03/903]T plate at 35Hz is 3.8g. It is therefore concluded that the analytical
model of Eq. (18) is capable to capture the characteristics of the chaotic snap-through
and continuous snap through events.

25


0 (t ) = 22 × t × 0.1g × sin(22 × 2π × t )
w
50 single-well intermittent continuous

Displacement of corner node (mm)

40 cross-well cross-well
30
20
10
0
-10
-20
-30
-40
-50
0.0 0.5 1.0 1.5 2.0 2.5
time (s)

(a) 100mm×100mm, [02/902]T plate

0 (t) = 35 × t × 0.1g × sin(35 × 2π × t )

w
30
single-well intermittent continuous cross-well
25
Displacement of corner node (mm)

cross-well
20
15
10
5
0
-5
-10
-15
-20
-25
-30
0.0 0.5 1.0 1.5 2.0 2.5
time (s)

(b) 100mm×100mm, [03/903]T plate

Fig. 18 The time domain responses of bi-stable plates excited by a sinusoidal input, predicted by
analytical model, i.e. Eq. (18). The excitation amplitude increases by 0.1g in each cycle, excitation


0 (t) = f × t × 0.1g × sin( f × 2π × t)
acceleration is w

The dynamic response of 100mm × 100mm, [02/902]T and [03/903]T plates at

different frequencies are predicted by both the analytical model given by Eq. (18), and
FEA. A constant damping coefficient of α = 100 is adopted. Sinusoidal excitations
with linearly increased amplitude are applied, that is w

0 (t) = 2t ×9.8×sin( f × 2π × t) .
Critical excitation amplitudes to induce cross-well vibration are determined according
to the vibration curve of the plate corner, as illustrated in Fig. 18.
The critical excitation amplitudes at different excitation frequencies are predicted
by both the analytical model and the FEA, and are compared in Fig. 19. For both
plates, the analytical and numerical methods predict identical varying rate of critical

26
excitation acceleration, and the lowest value. The overall curves given by Eq. (18) are
slightly shifted to the right side compared with the FEA curves. The analytical results
presented in Fig. 19 indicate that Eq. (18) slightly overestimates the stiffness of the
plates. Nevertheless, the results demonstrates that the analytical model is able to
accurately predict the nonlinearity of cross-ply bi-stable plates.

[02/902]T plate [03/903]T plate

12
7 Cross-well-analytical Cross-well-analytical
Cross-well-FEA Cross-well-FEA
10
Critical excitation amplitude (g)

Critical excitation amplitude (g)

6

5 8

4 6

3 4
2
2
1
0
18 20 22 24 26 28 30 25 30 35 40 45 50
Excitation frequency (Hz) Excitation frequency (Hz)

Fig. 19 Critical excitation acceleration needed to induce the cross-well vibrations for the 100mm
×100mm, [02/902 ]T and [03/903]T plates

8.5 Degree-of-Freedom of the analytical model and its efficiency

In the proposed analytical model, there are 17 terms in the assumed forms of w ( t ) ,
ε xm ( t ) and ε m
y ( t ) , as presented in Eqs. (9) and (10). The sixth order out-of-plane
displacement w ( t ) is assumed first. In order to test the minimum terms that are
required to obtain accurate results, the forms of ε xm ( t ) and ε m
y ( t ) are initially
expanded into complete polynomials. As it is aimed to minimize the unknowns in the
analytical model, the terms in ε xm ( t ) and ε ym ( t ) are deleted step by step in the testing
simulation. After many attempts, it was found that each term left in Eq. (10) is critical
and should be remained for obtaining accurate results. It was also interesting to note
that, in Eq. (10), ε xm ( t ) has identical terms of y as w ( t ) , and ε m
y ( t ) has identical
terms of x as w ( t ) .

For the illustration purpose, comparison of static load-displacement curves of the

100mm×100mm, [03/903]T plate predicted using different number of terms in the
analytical model. The numerical results are presented in Fig. 20. Concentrated loads
are applied on the corners of the plate, which is mounted at the central point. For the
15 degree-of-freedom analytical model, the c4 ( t ) x 2 y 2 and d 4 ( t ) x 2 y 2 terms in
ε xm ( t ) and ε m
y ( t ) , i.e. Eq. (10), are deleted. The c3 ( t ) y and d3 ( t ) x terms are
6 6

27
further deleted to reduce the degree-of-freedoms to 13. The 13, 15 and 17
degree-of-freedom in the analytical models have similar precision if the applied loads
are smaller or larger than the snap-through load. However, compared to the FEA
results, the prediction errors given by the analytical model on the snap through load is
increased when the degree-of-freedom is reduced. As the snap through phenomenon
occurs in cross-well vibration, inaccurate prediction of snap through will lead to
inaccurate dynamic analysis for bistable plates. In contrast, the 17 degree-of-freedom
analytical model predicts accurate snap through phenomenon of bistable plates, as it is
also validated by the FEA results. This numerical simulation further approves that
using 17 terms is the minimum degree-freedom for the analytical model proposed in
this work.

5
Out of plane displacement at plate corner (mm)

FEA
4 17-freedom
3 15-freedom
13-freedom
2
1 2N
2N

0 2N
2N

-1
2N
2N
-2
2N
-3 2N

-4
-5
-6
-7
0.0 0.5 1.0 1.5 2.0
Load (N)

Fig. 20 Static load-displacement curves of the 100mm×100mm, [03/903]T plate

In this study, the “NDSolve” function in Mathematica® is applied to solve the

nonlinear differential equations, i.e. Eqs. (18) and (22). The previous time domain
results presented in Figs. 7-18, each of which only requires less than one minute to
obtain the results using the proposed analytical model. In contrast, the time-domain
dynamic analysis using finite element method is very time-consuming. To obtain
more accurate results, the time step in the “Implicit, Dynamic” step should be
sufficiently small. In this study, the time step is set to be 0.0025s, therefore a set of
numerical results in time domain of 1 second needs 2000 iterations. A cross-well
vibration lasting for 2 seconds usually takes 20 to 40 minutes, depending on the mesh
size. In Fig. 19, the prediction of critical excitation acceleration by FEA costs much
longer computation time than the analytical model. This comparison of computational
time demonstrates the high efficiency provided by the analytical model.
9. Conclusion

The dynamics of bi-stable plates are complicated due to the nonlinearity and the
snap-through phenomena. Although several dynamic models have been reported,

28
 either the dynamic snap-through load was overestimated [24] or the dynamic response
at single points were predicted [22, 23]. An analytical dynamic model with only 17
unknown terms is established based on Rayleigh’s method and Hamilton’s principle in
this study. This established model can accurately predict the nonlinear dynamic
response of bi-stable plates. Moreover, this analytical model is able to predict the
nonlinear dynamic response over the entire region of bi-stable plates, rather than just
few specific points. The influence from the inertial mass and the coupling between a
bi-stable plate and the electrical shaker is also considered in the dynamic model. The
free vibration, excited vibrations including single-well, intermittent cross-well and
continuous cross-well, are predicted and studied. The predicted results of the
proposed analytical model for different vibration patterns of bi-stable plates coincide
very well with FEA results and experimental results.
The proposed reduced order analytical model for the nonlinear dynamic analysis
of composite plates possesses several advantages. Firstly, this proposed analytical
model provides a very efficient and effective means to predict the nonlinear dynamics
for cross-ply composite bi-stable plates. Secondly, this analytical model is very simple
to implement, provided that the material properties and dimensions are provided.
Thirdly, the procedure to predict the nonlinear dynamics of bi-stable plates using this
analytical model is a relatively simple extension to the static analysis. With only 17
unknown terms, the analytical model is not sophisticated enough to perform very
complicated dynamical analysis for bi-stable plates. Nevertheless, an extension of this
work to a more general model for angle-ply bi-stable composite plates is
straightforward.

Acknowledgement

Zhangming Wu would like to acknowledge the financial support from China’s 1000
young scholarship plan. Hao Li sincerely acknowledge the support of National
Natural Science Foundation of China Youth Fund, Grant No. 51605299

Appendix

2 8 1
u 0 (t ) = ( c(t ) + A−1 N sT1 ) x − a(t ) 2 x3 − a(t )a1 (t ) x5 − ( 8a1 (t )2 + 12a(t )a2 (t ) ) x 7
3 5 7
8 18
− a1 (t )a2 (t ) x9 − a2 (t )2 x11 + c1 (t ) xy 2 + c2 (t ) xy 4 + c3 (t ) xy 6
3 11
1 2 8 12
+ ( c4 (t ) − 4a(t )e(t ) ) x3 y 2 − e(t )2 x3 y 4 − a1 (t )e(t ) x5 y 2 − a2 (t )e(t ) x 7 y 2 (A.01)
3 3 5 7
 2 
− 2 ( a(t ) B*11 + b(t ) B*12 ) x −  4a1B*11 + e(t ) B*12  x3 − 6a2 (t ) B*11 x5
 3 
− (12b1 (t ) B*12 + 2e(t ) B*11 ) xy 2 − 30b2 (t ) B*12 xy 4

29
 2 8 1
v 0 (t ) = ( d (t ) + A−1 N sT1 ) y − b(t )2 y 3 − b(t )b1 (t ) y 5 − ( 8b1 (t )2 + 12b(t )b2 (t ) ) y 7
3 5 7
8 18
− b1 (t )b2 (t ) y 9 − b2 (t )2 y11 + d1 (t ) x 2 y + d2 (t ) x 4 y + d3 (t ) x6 y
3 11
1 2 8 12
+ ( d 4 (t ) − 4b(t )e(t ) ) x 2 y 3 − e(t )2 x 4 y 3 − b1 (t )e(t ) x 2 y 5 − b2 (t )e(t ) x 2 y 7 (A.02)
3 3 5 7
 2 
+ 2 ( b(t ) B*11 + a(t ) B*12 ) y +  4b1B*11 + e(t ) B*12  y 3 − 6b2 (t ) B*11 y 5
 3 
+ (12a1 (t ) B*12 + 2e(t ) B*11 ) x 2 y + 30a2 (t ) B*12 x 4 y

2 4 
γ xy0 = ( 4a(t )b(t ) + 2c1 (t ) + 2d1 (t ) ) xy +  c4 (t ) + 4d 2 (t ) + 8a1 (t )b(t ) + a(t )e(t )  x 3 y
3 3 
2 4   4 
+  d 4 (t ) + 4c2 (t ) + 8a(t )b1 (t ) + b(t )e(t )  xy 3 + 16a1 (t )b1 (t ) − e(t )2  x 3 y 3
3 3   3 
 24   24 
+  6d3 (t ) + 12a2 (t )b(t ) + a1 (t )e(t )  x5 y +  6c3 (t ) + 12a(t )b2 (t ) + b1 (t )e(t )  xy 5
 5   5  (A.03)
60
+24a2 (t )b1 (t ) x5 y 3 + 24a1 (t )b2 (t ) x3 y 5 + 36a(t ) 2 b(t ) 2 x5 y 5 + b2 (t )e(t ) xy 7
7
60
+ a2 (t )e(t ) x7 y + ( 24a1 (t ) B12* − 24b1 (t ) A−1B12* ) xy + 120a2 (t ) B*12 x3 y
7
−120b2 (t ) B*12 xy 3
Lx Ly
 ∂w(t ) ∂w(t ) 
M1 = ρ h ∫ 2
Lx∫ 2
Ly  dydx, i, j = 1 ~ 17 (A.04)
−
2
−
2  ∂X 1 (i ) ∂X 1 ( j ) 

 ∂(Π(t ) plate − WF (t )) 
K1 ( X1 ) =   , i = 1 ~ 17 (A.05)
 ∂X 1 (i) 
Lx Ly
 ∂w(t ) 
G = −ρh∫ 2
Lx ∫2
Ly  dydx, i = 1 ~ 17 (A.06)
−
2
−
2  ∂X 1 (i) 

 ∂w(t ) 
 1
Lx Ly 
m 0 2 
∂X 2 (i )
dxdy
 + ρ h ∫ Lx ∫ Ly  i, j = 2 ~ 18 (A.07)
2
M2 = 
 0 0 − − ∂w(t ) ∂w(t ) ∂w(t ) 
2 
∂X 2 (i) ∂X 2 ( j ) 
2
 ∂X 2 ( j )

  Lx Ly   Lx Ly  
 ∂w  2 , 2 , t  ∂w  2 , 2 , t  
M inertial1 = 4minertial       , i, j = 1 ~ 17 (A.08)
 ∂X 1 [i ] ∂X 1 [ j ] 
 
 

30
   Lx Ly  
 ∂w  , , t  
 1  2 2  
 ∂X 2 (i ) 
M inertial 2 = 4minertial   , i, j = 2 ~ 18 (A.09)
 ∂w  Lx , Ly , t   Lx Ly   Lx Ly  
∂w  , , t  ∂w  , , t 
  2 2   2 2   2 2 
 ∂X 2 ( j ) ∂X 2 (i ) ∂X 2 ( j ) 
 

  Lx Ly  
∂w , ,t
d (w0 (t ))   2 2  
2
F ( X 1 )inertial = −4minertial   , i = 1 ~ 17 (A.10)
dt 2  ∂X (i) 
 

Reference

1. Emam, S.A. and D.J. Inman, A Review on Bistable Composite Laminates for
Morphing and Energy Harvesting. Applied Mechanics Reviews, 2015. 67(6):
p. 060803.
2. Hyer, M.W., Some observations on the cured shape of thin unsymmetric
laminates. Journal of Composite Materials, 1981. 15(2): p. 175-194.
3. Hyer, M.W., The room-temperature shapes of four-layer unsymmetric
cross-ply laminates. Journal of Composite Materials, 1982. 16(4): p. 318.
4. Ren, L., A. Parvizi-Majidi, and Z. Li, Cured Shape of Cross-Ply Composite
Thin Shells. Journal of Composite Materials, 2003. 37(20): p. 1801-1820.
5. Dano, M.-L., M. Jean-St-Laurent, and A. Fecteau, Morphing of bistable
composite laminates using distributed piezoelectric actuators. Smart Materials
Research, 2012. 2012.
6. Daynes, S., P. Weaver, and J. Trevarthen, A morphing composite air inlet with
multiple stable shapes. Journal of Intelligent Material Systems and Structures,
2011: p. 1045389X11399943.
7. Pirrera, A., D. Avitabile, and P. Weaver, On the thermally induced bistability of
composite cylindrical shells for morphing structures. International journal of
solids and structures, 2012. 49(5): p. 685-700.
8. Kim, S.-W., et al. Towards a bio-mimetic flytrap robot based on a
snap-through mechanism. in Biomedical Robotics and Biomechatronics
(BioRob), 2010 3rd IEEE RAS and EMBS International Conference on. 2010.
IEEE.
9. Daynes, S., et al., On a Bistable Flap for an Airfoil, in 50th
AIAA/ASME/ASME/ASCE/AHS/ASC Structure, Structure Dynamics, and
Materials Conference. 2009, American Institute of Aeronautics and
Astronautics.
10. G.Diaconu, C., P. M.Weaver, and F. Mattioni, Concepts for morphing airfoil
sections using bi-stable laminated composite structures. Thin Walled

31
 Structures, 2007. 11: p. 1-13.
11. Arrieta, A.F., et al., A cantilevered piezoelectric bi-stable composite concept
for broadband energy harvesting, in SPIE Smart Structures and Materials+
Nondestructive Evaluation and Health Monitoring. 2013, International Society
for Optics and Photonics: San Diego, California p. 86880G-86880G-9.
12. Arrieta, A.F., et al., Broadband vibration energy harvesting based on
cantilevered piezoelectric bi-stable composites. Applied Physics Letters, 2013.
102(17): p. 173904.
13. Li, H., F. Dai, and S. Du, Broadband energy harvesting by exploiting
nonlinear oscillations around the second vibration mode of a rectangular
piezoelectric bistable laminate. Smart materials and structures, 2015. 24(4): p.
045024.
14. Arrieta, A.F., D.J. Wagg, and S.A. Neild, Dynamic Snap-through for Morphing
of Bi-stable Composite Plates. Journal of Intelligent Material Systems and
Structures, 2011. 22(2): p. 103-112.
15. Onur, B., et al., Dynamic control of a bistable wing under aerodynamic
loading. Smart Materials and Structures, 2013. 22(2): p. 025020.
16. Simsek, M.R. and O. Bilgen, Hybrid Position Feedback Controller for
Inducing Cross-Well Motion of Bistable Structures. AIAA Journal, 2016.
54(12): p. 4011-4021.
17. Zarepoor, M. and O. Bilgen, Constrained-Energy Cross-Well Actuation of
Bistable Structures. AIAA Journal, 2016. 54(9): p. 2905-2908.
18. Arrieta, A., et al., A piezoelectric bistable plate for nonlinear broadband
energy harvesting. Applied Physics Letters, 2010. 97(10): p. 104102.
19. Pellegrini, S.P., et al., Bistable vibration energy harvesters: A review. Journal
of Intelligent Material Systems and Structures, 2012. 24(11): p. 1303-1312.
20. Betts, D.N., et al., Nonlinear dynamics of a bistable piezoelectric-composite
energy harvester for broadband application. The European Physical Journal
Special Topics, 2013. 222(7): p. 1553-1562.
21. Harne, R.L. and K.W. Wang, A review of the recent research on vibration
energy harvesting via bistable systems. Smart Materials and Structures, 2013.
22(2): p. 023001.
22. Arrieta, A., S. Neild, and D. Wagg, Nonlinear dynamic response and modeling
of a bi-stable composite plate for applications to adaptive structures.
Nonlinear Dynamics, 2009. 58(1): p. 259-272.
23. Arrieta, A.F., S.A. Neild, and D.J. Wagg, On the cross-well dynamics of a
bi-stable composite plate. Journal of Sound and Vibration, 2011. 330(14): p.
3424-3441.
24. Diaconu, C.G., P.M. Weaver, and A.F. Arrieta, Dynamic analysis of bi-stable
composite plates. Journal of Sound and Vibration, 2009. 322(4): p. 987-1004.
25. Vogl, G.A. and M.W. Hyer, Natural vibration of unsymmetric cross-ply
laminates. Journal of Sound and Vibration, 2011. 330(20): p. 4764-4779.

32
26. Arrieta, A.F., et al., Low order model for the dynamics of bi-stable composite
plates. Journal of Intelligent Material Systems and Structures, 2011. 22(17): p.
2025-2043.
27. Taki, M.S., et al., Dynamic responses of cross-ply bi-stable composite
laminates with piezoelectric layers. Archive of Applied Mechanics, 2016.
86(6): p. 1003-1018.
28. Peeters, L., P. Powell, and L. Warnet, Thermally-induced shapes of
unsymmetric laminates. Journal of Composite Materials, 1996. 30(5): p.
603-626.
29. M.W.Hyer, Calculations of the Room-Temperature Shapes of Unsymmetric
Laminates. Composite Materials, 1981. 15: p. 296-310.
30. Aimmanee, S. and M. Hyer, Analysis of the manufactured shape of
rectangular THUNDER-type actuators. Smart materials and structures, 2004.
13: p. 1389.
31. Pirrera, A., D. Avitabile, and P. Weaver, Bistable plates for morphing
structures: A refined analytical approach with high-order polynomials.
International journal of solids and structures, 2010. 47(25-26): p. 3412-3425.
32. Ashton, J.E., Approximate Solutions for Unsymmetrically Laminated Plates.
Journal of Composite Materials, 1969. 3: p. 189-191.
33. Gigliotti, M., et al., Predicting loss of bifurcation behaviour of 0/90
unsymmetric composite plates subjected to environmental loads. Composite
Structures, 2012.
34. Arrieta, A., D. Wagg, and S. Neild, Dynamic snap-through for morphing of
bi-stable composite plates. Journal of Intelligent Material Systems and
Structures, 2011. 22: p. 103-112.

33