1. science
2. mathematics
3. algebra

Proceedings of the 9th International Conference on Structural Dynamics, EURODYN 2014
Porto, Portugal, 30 June - 2 July 2014
A. Cunha, E. Caetano, P. Ribeiro, G. Müller (eds.)
ISSN: 2311-9020; ISBN: 978-972-752-165-4
Reduced-order models for the analysis of non-linear vibrations
of variable stiffness composite laminated plates
Hamed Akhavan, Pedro Ribeiro
Department of Mechanical Eng., Faculty of Engineering, University of Porto, R. Dr. Roberto Frias, 4200-465 Porto, Portugal
email: hamed@fe.up.pt, pmleal@fe.up.pt
ABSTRACT: This study focuses on reduced-order models to analyse non-linear oscillations of rectangular variable stiffness
composite laminated (VSCL) plates. Each lamina of these VSCL plates has curvilinear fibres that are distributed by shifting
along one of the Cartesian directions. Two model reduction approaches are applied: modal reduction alone and static
condensation followed by modal reduction. In both cases, a p-version finite element that takes into account the effect of large
amplitude displacements in the framework of third-order shear deformation theory (TSDT) is used to compute the linear modes,
which are used as a basis. In the first approach, a reduced set that contains bending and in-plane modes is selected. The mode
selection is discussed and it is confirmed that, also in this TSDT based model and in the non-linear regime, in-plane modes are
essential for accuracy. In the second approach, the in-plane inertia is neglected in the full model and static condensation is
applied, so that the influence of in-plane displacements is kept; only then is modal reduction implemented. This model has the
advantage of leading to a smaller number of degrees of freedom and of automatically taking care of the in-plane displacements,
but loses the effect of in-plane inertia. The accuracy and computational performance of both models in the definition of periodic
responses of VSCL plates is compared. For that purpose, the laminates are subjected to transverse harmonic loads and the
frequency responses are determined using the shooting method, with fifth-order Cash-Karp Runge-Kutta method and adaptive
stepsize control. For excitation frequencies around the first natural frequency, the amplitudes of vibration, the time histories and
the stability status of the diverse solutions are compared.
KEY WORDS: Variable stiffness; Curvilinear fibres; Non-linear; Periodic vibration; Reduced-order model.
1
INTRODUCTION
Variable stiffness composite laminated (VSCL) plates are a
general type of composite plates, in comparison to constant
stiffness composite laminated (CSCL) plates, see Ref. [1]. In
the VSCL plates investigated here, fibres follow a curvilinear
path (their path is straight in CSCL plates). Mechanical
properties of composite plates with curvilinear fibres can
change with fibre paths. Using curvilinear fibres instead of
straight ones can improve these properties in specific
structures (e.g. at the corners of aircraft fuselage windows).
Some favourable factors of VSCL against CSCL plates in
static non-linear analyses and free linear vibration are
mentioned, for example, in Refs. [2-5]. Few works [6-8]
addressed free and forced non-linear vibration of VSCL plates
by using thin plate theory (CPT) and first-order shear
deformation theory (FSDT). In references [2-8], the p-version
finite element, with hierarchic shape functions, is used. In
references [6, 7] the harmonic balance method (HBM) is
applied, whereas in ref. [8] Newmark method is used. A
comprehensive review on VSCL plates and methods used to
analyse them can be found in Ref. [1].
In the current paper, we intend to apply a p-version finite
element model, based on a third-order shear deformation
theory (TSDT), to VSCL plates with shifted curvilinear fibres
[9]. The finite element model may be condensed statically
[10] and is reduced using modal coordinates [11]. The choice
of modes is analysed in two model reduction procedures.
Periodic solutions of the plate under harmonic loading, in the
non-linear regime, are calculated using the shooting method.
The monodromy matrix, which is a by-product of the shooting
method, is employed to define the stability of the solutions.
The response of VSCL plates to harmonic excitations with
frequencies close to their fundamental frequencies is
investigated. For that purpose, the frequency response curves
are given, as well as phase plots and time histories. Fourier
spectra are used to evaluate the harmonic content of the
solution.
2
2.1
PERIODIC VIBRATION OF VSCL PLATES
p-Version finite element model
Rectangular, symmetric about the middle plane, clamped
laminates are displayed in Figure 1, as well as Cartesian
coordinates with origin at the plate centroid. Fibres in the ith
lamina can be either straight ( Τ 0 = Τ 1 ), where Τ 0 and Τ 1 are
the fibre angles at the centre and vertical edges, or curvilinear
( Τ 0 ≠ Τ 1 ). The fibre angles in the ith lamina change as
θi  x   2 Τ1 - Τ 0  x a  Τ 0 . Then, a configuration like [〈 Τ 0 ,
Τ 1 〉1,…, 〈 Τ 0 , Τ 1 〉i,…, 〈 Τ 0 , Τ 1 〉k]sym represents a 2k-layered
symmetric laminate with different fibre angles in each layer.
Displacement field, ix, y, z, t , i  u, v, w , following a
third-order shear deformation theory (TSDT) in x, y, z
coordinates, respectively, is (with
)


u  u 0  z x  cz 3  x  w0 x ,
1803
Proceedings of the 9th International Conference on Structural Dynamics, EURODYN 2014


v  v 0  z y  cz 3  y  w0 y ,
(1)
w  w0 ,
in which i x, y, t , i  u, v, w are displacements at the midplane and i x, y, t , i  x, y are rotations around axis y and x.
0
a
T0
T1
h
z
x
z
Figure 1. Up is a ply with a reference curvilinear fibre path
and down is a four-layer symmetric laminate. T0 and T1 have
negative values in this picture.
A p-version finite element method (FEM) with hierarchical
basis functions that can be applied to thin and thick plates,
because it is based on a TSDT [2-4], is used. The
displacement components are related to the generalized
displacements by
where
0
T
Nu
0
0
0
0
0
T
Nw
0
0
N i x, y , i  u, w,  x ,  y
0
0
0
N
x T
0
  q u t  


  q v t  

  q w t 





q
t
  
 T  q t 
N 
  

0
0
0
0

 yz  1  3cz  y  w 0 y ,
(3)
2
 xz  1  3cz 2  x  w 0 x .
0
0
Q44
0
0
0
0
0
Q55
0
0  1 
0    2 
 
0   23 

0   13 
 
Q66   12 
(4)
where subscripts 1, 2, and 3 represent, respectively, the fibre
orientation, the in-plane orientation perpendicular to 1, and
the transverse direction. Stresses and strains in plate
coordinates x and y are found using transformation matrices
(see Ref. [14]). Comparing constant stiffness to variable
stiffness composites, it is worthy to say that due to curvilinear
fibres in VSCL, directions 1 and 2 are not anymore constant.
Defining virtual works of inertia, internal and external
forces and then applying them in the principle of virtual work,
leads to the equations of motion, as shown below (details on
the matrices are given in [3]).
y
 u 0  N u T
 0 
v   0
 0 
w   0
  
 x  0


 y 
  0

 1  Q11 Q12
  Q
 2   21 Q22
0
 23    0
   0
0
 13  
 12   0
0
x
y

 w 0 x w 0 y  2cz 3  2 w 0 xy ,
Each ply with curvilinear fibre is a non-homogenous
orthotropic fibrous composite layer, in which stresses and
strains are related as [14]
T1
b
 xy  u 0 y  v 0 x  z  cz 3  x y   y x 
(2)
x
y
y
are vectors composed of in-
plane, out-of-plane and rotational shape functions, found for
instance in [12, 13]. The strain-displacement relations can be
approximated using the von Kármán strains by
 x  u 0 x  z  cz 3  x x  cz 3  2 w0 x 2  w0 x  2 ,
2
 y  v 0 y  z  cz 3  y y  cz 3  2 w0 y 2  w0 y  2 ,
M 11
0

M 22




 sym

 K 11
K 12
L
L

22
K
L




 sym

0
0
M 33
0
0
M 34
M 44
 u t  
0  q

 q

0    v t  
  t 
M 35   q
w


  t 
0  q

  t 
M 55  q



q
t



0
u
  q t  
0  v 
  q t 
K 35
L  w

K L45  q  t 
 q t 
K 55
L  

x
y
0
0
K 33
L
0
0
K 34
L
K L44
x
y

0
0
K 13
NL q w t 

23
0
0
K
NL q w t 

31
32
33

 K NL q w t  K NL q w t  K NL q w t 

0
0
0


0
0
0

 0 
 0 


 f w t .
 0 


 0 
0
0
0
0
0
0  q u t  


0  q v t  


0  q w t 

0 q  t 
0 q  t 
x
y
(5)
2
 z  0,
1804
The vector on the right-hand side of equation (5) represents
generalized external forces.
For example, taking 7 shape functions in x and y directions
– equal to a vector of forty nine 2-D shape functions – for
each of five variables in equation (2), leads a full model in
Proceedings of the 9th International Conference on Structural Dynamics, EURODYN 2014
equation (5) with n=245 DOF (hereafter, n represents the
number of degrees of freedom, or “order”, of the full model).
2.2
Static condensation
Due to the clamped boundaries of the plate, absence of inplane external force and the fact that the transverse deflection
is not bigger than two times the plate’s thickness, the in-plane
displacements are very small and the corresponding in-plane
inertia is not vital to the solution. Then the in-plane inertia can
be eliminated from the formulation. This elimination process
is known as static condensation [10] and reduces the order of
the problem. Although it generally changes the bandwidth of
the stiffness matrix, in the present TSDT formulation – and as
occurs with Kirchhoff [15] or FSDT [12] approaches – the
bandwidth of the linear matrices are not increased by static
condensation. The equations of motion (5) are condensed as
 w t  K 33
M 35   q
K 34
L
L
 

  t   
0  q
K L44
 t   sym
M 55  q
   
K NL q w t  0 0  q w t  f w t 

 

 
0 0 q  t    0 .

sym
0 q  t   0 


M 33


 sym

M 34
M 44
x
y
with
T

Shooting method
To apply the shooting method [16], the system of m secondorder ordinary differential equations of motion (11) is
transformed to a system of 2m first-order ordinary differential
equations using q m  y ,
0

M
3
1
K  K
T K
K 23
 
NL 
sym K  K
11
L
12
L
22
L
13
NL
23
NL

. (7)

If the full model is of n=245 DOF, the condensed model
would be of n=147 DOF. Both equation (5) and (6) may
simply be written as
t   K Lqt   K NL qt qt   f t .
Mq
2.3
2.4
  y  0 
0
M   y  - M
     .
    
0  q m   0 K L  K NL  q m  p
(12)
stiffness terms. Apart from the improvements in the
integration of first-order differential equations (12), the
procedure here applied is similar to the one described in Refs.
[17, 18].
x

non-linear stiffness matrix K NL is not. In the above example,
the full model was of order n=245 DOF, taking m=10 modes
results in the reduced modal model of m=10 DOF.
Shooting with fifth-order Runge-Kutta method with stepsize
adjustment algorithm based on the embedded fourth-order
  q w t 
K 35
Runge-Kutta formulas [16] (with constant parameters found
L


by Cash and Karp) is applied, starting by using as first guess
K L45  q x t 
55  

y 
K L  q y t 
of   the solution of equation (12) only with linear


(6)
q m 
y
13
K NL  K 33
NL  2 K NL
linear stiffness matrix K L are diagonal, but the generalized
(8)
Analysis in modal coordinates
If we know that the excitation of the VSCL plate centres
around a specific frequency, many of the modes probably will
not be excited and we can assume the forced vibration to be
the superposition of only a few modes [11], let’s say m modes.
The truncated modal matrix  nm assembled by m normal
modes ( i as eigenvectors of linear matrix) can appear as
  i ... m , m  n
(9)
Then the displacements can be written as
qt   q m t 
(10)
Pre-multiplying equation (8) by the transpose  T reduces the
n-DOF system to one with a dimension equal to the number of
modes used, m-DOF. Subsequently, instead of solving n
coupled equations (8), we need to solve the m coupled
equations represented by
 m t   K L q m t   K NL q m t q m t   pt .
Mq
(11)
We are using the modes of the linear system in a non-linear
problem; therefore, the generalized mass matrix M and the
3.1
NUMERICAL TESTS AND VALIDATION
How to select modes in two reduced-order models
In the first approach, reduced-order model I (model I, in
brief), the full model of equation (5) is transferred to a
reduced set of modal coordinates. The modes to include in the
truncated modal matrix  , should be the more influential
bending modes along with influential in-plane modes. Here,
expression “bending mode” refers to normal modes of the
linear system that actually contain transverse and rotational
generalized displacements; on the other hand, “in-plane
mode” identifies the normal modes with only membrane
generalized displacements. It is obvious from the mass and
linear stiffness matrices of equation (5), of model I, that
membrane and bending modes are uncoupled. We note that
shear is also present in the model, and we may have bending
coupled with shear, and modes where shear is more important;
for the sake of simplicity, we will designate modes where
bending is important as “bending modes”.
In the second approach, reduced-order model II (model II),
the condensed model of equation (6), which has not in-plane
inertia but is still affected by the influence of in-plane
displacements, is transformed to modal coordinates. In this
second approach, only bending modes are selected.
The procedure of selecting influential modes is programmed
and implemented as an in-house code in FORTRAN. In
reduced-order model I, firstly all modes of vibration of the
linear problem, including in-plane normal modes, are found.
Then, equation (11) is solved in modal coordinates by the
shooting method, using only one bending mode. A symmetric
distributed transverse load with a fixed excitation frequency is
applied, because this is the type of load considered in this
work. The process is repeated for all non-in-plane modes, one
by one. As an example, taking a vector of forty nine 2-D
1805
Proceedings of the 9th International Conference on Structural Dynamics, EURODYN 2014
shape functions for each variable of equation (2), this step
(Step 1) repeats for each of 3×49 modes, where transverse
displacements are coupled with cross section rotations. In this
step, those bending modes that result in bigger central
deflection are candidates to enter the next step.
In the next step, Step 2, equation (11) is again solved, but
with two or three modes. One of these modes is the bending
mode selected from the previous step, while the other one or
two modes is/are selected among the in-plane modes of the
linear problem. This stage may be repeated for all candidate
bending mode selected in Step 1, combined with different inplane modes. Again, like in the end of step one, groups of
bending and in-plane modes with bigger deflections are
selected. The procedure ends when enough influential modes
Table 1. Frequency ratio
±0.2
±0.6
±1
3.2
are found. Here, the word enough means a number of modes
that gives convergence in deflection computation.
To select normal modes for reduced-order model II, the
same procedure is applied, but the first step is sufficient, due
to the lack of in-plane modes in linear model II. In this model,
the influence of in-plane displacements is considered in the
reduced non-linear matrix.
Then, the analysis to find forced responses of vibration of
VSCL plates is performed by applying enough influential
modes, which were selected by one of the aforementioned
procedures, resulting in two types of reduced-order models. In
the following section the two types of reduced order
procedures are compared.
against amplitude of vibration for an isotropic plate.
HFEM and
shooting
Ref. [18]
FEM
without IDI
Ref. [19]
Elliptic
function
[20, 21]
Newmark
method
Model I
(8 modal
coordinates)
Model I
(15 modal
coordinates)
Model II
(3 modal
coordinates)
Model II
(15 modal
coordinates)
0.2707
1.4387
0.8950
1.2114
1.0777
-
0.1180
1.4195
0.8905
1.2083
1.0700
1.2429
0.1200
1.4195
0.8951
1.2117
1.0822
1.2540
0.2400
1.4394
0.8961
1.0825
-
0.131
1.418
0.914
1.225
1.126
1.298
0.2797
1.4418
0.9044
1.2674
1.1007
1.3332
0.1338
1.4202
0.8905
1.2106
1.0744
1.2494
0.1203
1.4189
0.8904
1.2102
1.0748
1.2498
Comparison of two reduced-order models
First, a square isotropic plate with fully clamped boundaries is
analysed. Table 1 compares the amplitude ratios
of
transverse central deflection of an undamped steel plate, as a
function of excitation frequency, calculated in Refs. [18-21] –
where models and methods of solution different from the ones
here proposed were applied – with: the proposed reducedorder models I and II and the shooting method; the full model
and Newmark method. Newmark method [22] is applied to
the full model, equation (5), to find deflections without
transformation to the modal coordinates, providing reference
solutions of the present model. Two sets of results are taken
from Ref. [19]. In one of them, an approach close to a
harmonic balance procedure with one harmonic is applied to
an h-version finite element model. In the other, following [20]
and [21], a one mode approach (one linear mode, from a
model with in-plane deformation but without inertia) is
employed to obtain a Duffing type equation that is exactly
solved using elliptic functions, while in Ref. [18] the shooting
method was applied. The plate is square with length a=500
mm and thickness h=2.0833 mm, and is subjected to a
harmonic distributed load
, ,
N/m2.
In Ref. [18], the p-version finite element method (HFEM) has
just three out-of-plane and six in-plane unidimensional shape
functions. In the table, IDI stands for in-plane displacement
and inertia. In the full model plus Newmark method, 245 DOF
are used; the same 245 DOF full model is employed to obtain
the reduced models that lead to the four columns with results
from the present paper, solved with the shooting method.
Increasing the number of modal coordinates from 3 to 15 in
model II, the frequencies change slightly. With either 3 or 15
modal coordinates, model II frequencies are close to the ones
provided by the full model, Newmark solution, with the
exception of the far away from resonance case Wmax/h=0.2.
Also, convergence to the full model solutions was not yet
achieved in model I using 15 modal coordinates; the values
are generally close but marginally worse than the ones of
model II. It should be mentioned that the run-time of model I,
with the same number of modal coordinates, is larger than the
run-time of model II.
Table 2. Comparison of the non-linear frequency ratio
in function of the amplitude of vibration
angle-ply plate, by the full model plus Newmark method.
0.2016
0.5997
1.0005
1806
p-version FEM CPT
[24]
0.7340
1.0085
1.1756
0.2
0.6
1
p-version FEM
FSDT [23]
0.7268
1.005
1.171
h-version FEM
CPT [25]
0.7219
1.0085
1.1787
, for a CSCL
Newmark method
245 DOF
0.7275
1.0097
1.1814
Proceedings of the 9th International Conference on Structural Dynamics, EURODYN 2014
In Table 2, a comparison is carried out with data available
on Refs. [23] and [24] (computed using the p-version FEM,
FSDT and thin plate theory, or CPT, theories, respectively)
and with data computed by the h-version FEM using CPT in
Ref. [25]. The full model (245 DOF) is solved in the timedomain by Newmark method [22], with a small stiffnessproportional damping parameter
, [26]. The CSCL
plate analysed is an angle-ply three-layer plate [45, -45, 45],
with the following mechanical and geometric characteristics:
,
,
,
,
,
,
The amplitude of deflection at centre of
the plate with a sinusoidal wave with 3454.665
is
given. The comparison of Table 2 allows us to verify the full
TSDT model.
Table 3. Comparison of the amplitude of vibration
of a CSCL angle-ply plate, in function of frequency,
different combinations of modes in modal coordinates and a full model by Newmark method.
0.7275
1.0097
1.1814
by Newmark method,
245 DOF
0.2
0.6
1.0
Bending modes
20 modes
30 modes
0.1935
0.1935
0.5089
0.5091
0.8244
0.8239
by Reduced-order model I
Including bending/in-plane modes
40 modes
19 modes
0.1935
0.1979
0.5092
0.5925
0.8243
0.9706
Table 4. Influence of bending modes and bending/in-plane modes in the amplitude of vibration
in function of excitation frequency .
Frequency
of external
load,
Deflection
ratio by
Newmak
method, full
model,
600
0.1175
650
0.1275
700
0.1404
750
0.1572
800
0.1799
850
0.2115
900
0.2570
950
0.3233
1000
0.4160
1050
0.5315
1100
0.6572
1150
0.7839
1200
0.9062
1250
1300
-
, using
, for an angle-ply plate
, reduced-order model I
5 modes
0.1174
(-0.04%)
0.1273
(-0.15%)
0.1401
(-0.24%)
0.1566
(-0.40%)
0.1780
(-1.05%)
0.2074
(-1. 94%)
0.2516
(-2.09%)
0.3120
(-3.50%)
0.3929
(-5.55%)
0.4886
(-8.07%)
0.5952
(-9.44%)
0.7046
(-10.1%)
0.7920
(-12.6%)
0.9086
1.0041
Only bending modes are taken
10 modes
15 modes
20 modes
0.1172
0.1173
0.1173
(-0.22%)
(-0.19%)
(-0.16%)
0.1271
0.1272
0.1272
(-0.29%)
(-0.23%)
(-0.22%)
0.1398
0.1399
0.1399
(-0.46%)
(-0.38%)
(-0.36%)
0.1562
0.1564
0.1565
(-0.61%)
(-0.50%)
(-0.46%)
0.1783
0.1786
0.1787
(-0.89%)
(-0.71%)
(-0.66%)
0.2087
0.2091
0.2094
(-1.33%)
(-1.11%)
(-1.01%)
0.2513
0.2524
0.2528
(-2.20%)
(-1.78%)
(-1.65%)
0.3112
0.3137
0.3143
(-3.75%)
(-2.97%)
(-2.77%)
0.3927
0.3967
0.3979
(-5.61%)
(-4.65%)
(-4.35%)
0.4909
0.4977
0.4998
(-7.64%)
(-6.36%)
(-5.96%)
0.5976
0.6069
0.6107
(-9.07%)
(-7.66%)
(-7.08%)
0.7054
0.7179
0.7224
(-10.0%)
(-8.42%)
(-7.85%)
0.8103
0.8258
0.8312
(-10.6%)
(-8.87%)
(-8.27%)
0.9087
0.9298
0.9363
1.0053
1.0301
1.0375
25 modes
0.1174
(-0.11%)
0.1273
(-0.15%)
0.1400
(-0.25%)
0.1567
(-0.32%)
0.1791
(-0.46%)
0.2100
(-0.70%)
0.2540
(-1.15%)
0.3172
(-1.88%)
0.4033
(-3.06%)
0.5097
(-4.11%)
0.6246
(-4.96%)
0.7406
(-5.53%)
0.8532
(-5.84%)
0.9622
1.0670
19
bending
and
in-plane
modes are
taken
0.1178
(0.23%)
0.1280
(0.37%)
0.1457
(3.76%)
0.1571
(-0.04%)
0.1798
(-0.05%)
0.2114
(-0.04%)
0.2569
(-0.03%)
0.3231
(-0.08%)
0.4163
(0.07%)
0.5320
(0.10%)
0.6575
(0.05%)
0.7839
(0.00%)
0.9065
(0.03%)
1.0241
1.1371
1807
Proceedings of the 9th International Conference on Structural Dynamics, EURODYN 2014
The full model, to some extent validated in Table 2, is used
as reference in Table 3, in another test of the model reduction
procedures. The plate is the one of Table 2. In Table 1,
increasing modal coordinates from 8 to 15 in model I, so that
more in-plane modes are considered, the relation between
frequency and amplitude of vibration was more accurately
computed. Table 3 also shows that, in reduced-order model I
and in the non-linear regime, in-plane modes are essential for
accuracy. Here, four analyses of non-linear frequencies of
undamped plate [45, -45, 45] are performed using modal
coordinates (of a model with 245 DOF): three of them include
only the first 20, 30, and 40 bending modes of vibration and
the last analysis involves a combination of 19 bending and inplane modes of vibration, selected by the technique written in
section 3.1.
In reduced-order model I, the results demonstrate that
increasing the number of bending modes alone cannot
guarantee the accuracy in the non-linear regime, but taking
less bending modes along with in-plane modes can profoundly
improve the precision of the model. Again, it is worthy to note
that reduced-order model II (that excludes in-plane inertia)
already considers the in-plane displacements, during static
condensation.
For the same plate and loading, a smaller model (using 9 bidimensional shape functions, i.e. 45 DOF) is applied in Table
4 to compare the amplitudes of vibration computed by
Newmark method (full model) with the ones computed by the
shooting method (reduced-order model I). Again, different
reduced-order models with different number of bending
normal modes in modal coordinates (5 to 25 modes) are
considered and compared with a model that incorporates a
combination of 19 bending/in-plane modes chosen as
explained in section 3.1. In this table, the numbers in
parenthesis represent the percentage of the difference with
respect to Newmark results (full model). The table verifies
that deflections using reduced-order model I and taking into
account only bending normal modes are always smaller than
reference deflections. Away from resonance, models only
with bending modes give accurate results. However, when the
non-linearity is more important, deflections with 19 bending
and in-plane modes are in a better agreement with the
reference data, than deflections computed with more bending
modes, but without in-plane modes.
a)
1,5
b)
CSCL - T1=45 (Stable)
VSCL - T1=45 (Unstable)
VSCL - T1=0 (Stable)
CSCL - T1=0 (Unstable)
VSCL - T1=25 (Stable)
VSCL - T1=25 (Unstable)
1
1,2
CSCL - T=45
VSCL - T=0
VSCL - T=25
W/h
0,8
0,4
W/h
0
0,5
-0,4
-0,8
0
-1,2
0
0,5
1
ω/ωl1
1,5
0
CSCL - T=45
𝑊
d)
0,8
VSCL - T=25
10
𝑊
/
0
0,5
𝑖/
0,6
0,4
0
-0,5
0,015
CSCL - T1=45
VSCL - T1=0
VSCL - T1=25
𝑊_
VSCL - T=0
-1
0,01
Time (s)
c)
20
0,005
1
0,2
-10
𝜔_𝑖∕
0
-20
0
1
2
3
4
5
6
7
8
9
10
Figure 2. Plate’s frequency response: a) maximum transverse displacement at centre of the plate with respect to frequency with
stable and unstable solutions; b), c) and d) displacements along one cycle, phase plot, and Fourier series at excitation frequency
equal to 0.378 , stable solution.
1808
Proceedings of the 9th International Conference on Structural Dynamics, EURODYN 2014
3.3
Forced vibration of VSCL plates with reduced-order
model
A test case of a clamped three-layered VSCL plate with fibre
configuration [〈45,T〉,〈-45,-T〉,〈45,T〉] is analysed, taking
constant Τ 0  45 in all layers and different Τ 1  T as 0, 25, 45
(while Τ 1  45 exemplifies the angle-ply CSCL plate [45,45,45]). The plate characteristics are:
,
,
,
,
,
,
The thickness of all layers is the
same. Reduced-order models with 6 modal coordinates are
used in the analysis; the corresponding full models had 245
DOF (49 bi-dimensional shape functions per displacement and
cross section rotation component) and the statically condensed
models had 147 DOF.
The responses of two VSCL plates, T=0 and 25, and a
CSCL plate, T=45, to a uniform distributed harmonic
excitation of amplitude 3454.665
are represented in
Figure 2. The plates are undamped. Instabilities when the
Floquet multiplier crosses the unit circle through +1 are found
and represented by a dashed line. The displacements along
one cycle, phase plots and Fourier spectra show that, at some
excitation frequencies, higher harmonics have a large
influence in the response. The width of the resonant region
shrinks slightly by changing T from 45 (i.e. CSCL plate) to 0
or 25 (i.e. VSCL plate). If the excitation frequency is equal to
0.378 , a similar pattern of displacements along a period of
vibration and similar phase plots are observed for CSCL and
VSCL plates, but the vibration period changes from plate to
plate.
4
CONCLUSION
A procedure based on the shooting method and on a Thirdorder Shear Deformation p-version finite element was
proposed to investigate the geometrically non-linear forced
vibration of variable stiffness composite plates. The method
can find instabilities and periodic oscillations with higher
harmonics, when they occur, for instance due to the
appearance of superharmonics in the response. The technique
takes advantage of two reduction methods, namely static
condensation and modal coordinates, to decrease the order of
matrices. If static condensation is not applied, then in-plane
modes must be considered in the reduced order model, to
achieve accurate results in the non-linear regime. Static
condensation permits to take in-plane displacements into
account indirectly, so that one only has to select bending
modes.
ACKNOWLEDGMENTS
This investigation is carried out under supports of Portuguese
Science and Technology Foundation (FCT) under PhD grant
SFRH/BD/81707/2011
and
project
PTDC/EMEPME/098967/2008, which are gratefully acknowledged.
[3] H. Akhavan, P. Ribeiro and M.F.S.F. de Moura, Large deflection and
stresses in variable stiffness composite laminates with curvilinear fibres,
International Journal of Mechanical Sciences, 73, 14-26, 2013.
[4] H. Akhavan, P. Ribeiro and M.F.S.F.de Moura, Composite laminates
with linearly varying fiber angles under static and dynamic loads,
Proceeding of 54th AIAA/ASME/ASCE/AHS/ASC Structures, Structural
Dynamics, and Materials Conference and co-located Events, Boston,
MA, 2013, AIAA 2013-1565.
[5] S. Yazdani, P. Ribeiro and J.D. Rodrigues, A p-version layerwise model
for large deflection of composite plates with curvilinear fibres,
Composite Structures, 108, 181-190, 2014.
[6] A. Houmat, Nonlinear free vibration of laminated composite rectangular
plates with curvilinear fibers, Composite Structures, 106, 211-224, 2013.
[7] P. Ribeiro, Non-linear free periodic vibrations of variable stiffness
composite laminated plates, Nonlinear Dynamics, 70, 2, 1535-1548,
2012.
[8] P. Ribeiro and H. Akhavan, Non-linear vibrations of variable stiffness
composite laminated plates, Composite Structures, 94, 8, 2424-2432,
2012.
[9] C. Waldhart, Analysis of tow-placed, variable-stiffness laminates, M.Sc.
thesis, Blacksburg, Virginia Tech, 1996.
[10] L. Meirovitch, Computational methods in structural dynamics, Sidthoff
& Noordhoff, Rockville, MD, US, 1990.
[11] W.T. Thomson and M.D. Dahleh, Eds., Theory of vibration with
applications, Prentice-Hall, Inc., New Jersey, US, fifth edition, 1998.
[12] P. Ribeiro, A hierarchical finite element for geometrically non-linear
vibration of thick plates, Meccanica, 38, 1, 117-132, 2003.
[13] P. Ribeiro and M. Petyt, Nonlinear vibration of plates by the hierarchical
finite element and continuation methods, International Journal of
Mechanical Sciences, 41, 4-5, 437–459, 1999.
[14] J.N., Reddy, Ed., Mechanics of laminated composite plates: theory and
analysis, CRC Press, Boca Raton, US, 1997.
[15] W. Han and M. Petyt, Linear vibration analysis of laminated rectangular
plates using the hierarchical finite element method– I. Free vibration
analysis, Computers & Structures, 61, 4, 705-712, 1996.
[16] W.H. Press, S.A. Teukolsky, B.P. Flannery and W.T. Vetterling, Eds.,
Numerical recipes in Fortran, Cambridge University Press, Cambridge,
New York, second edition, 1992.
[17] P. Ribeiro, Non-linear forced vibrations of thin/thick beams and plates by
the finite element and shooting methods, Computers & Structures, 82,
17-19, 1413-1423, 2004.
[18] P. Ribeiro, Periodic vibration of plates with large displacements, AIAA
Journal, 40, 1, 185-188, 2002.
[19] C. Mei and K. Decha-Umphai, Finite element method for nonlinear
forced vibrations of rectangular plates, AIAA Journal, 23, 7, 1104-1110,
1985.
[20] J.G. Eisley, Nonlinear vibration of beams and rectangular plates,
Zeitschrift fur angewandte Mathematik and Physik, 15, 167-175, 1964.
[21] C.S. Hsu, On the application of elliptic functions in nonlinear forced
oscillations, Quarterly on Applied Mathematics, 17, 393-407, 1960.
[22] K.J. Bathe and E.L. Wilson, Eds., Numerical methods in finite element
analysis, Prentice-Hall, cop., Englewood Cliffs, US, 1976.
[23] P. Ribeiro, Forced periodic vibrations of laminated composite plates by a
p-version, first order shear deformation, finite element, Composites
Science and Technology, 66, 11-12, 1844-1856, 2006.
[24] P. Ribeiro, and M. Petyt, Non-linear vibration of composite laminated
plates by the hierarchical finite element method, Composite Structures,
46, 3, 197-208, 1999.
[25] C.K. Chiang, C. Mei, and C.E. Gray, Finite element large-amplitude free
and forced vibrations of rectangular thin composite plates, Journal of
Vibration and Acoustics-Transactions of the ASME, 113, 3, 309-315,
1991.
[26] S. Adhikari, Damping models for structural vibration, Ph.D. thesis,
Engineering Department, Cambridge University, Cambridge, UK, 2000.
REFERENCES
[1] P. Ribeiro, H. Akhavan, A. Teter and J. Warminski, A review on the
mechanical behaviour of curvilinear fibre composite laminated panels,
Journal of Composite Materials, DOI: 10.1177/0021998313502066.
[2] H. Akhavan and P. Ribeiro, Natural modes of vibration of variable
stiffness composite laminates with curvilinear fibers, Composite
Structures, 93, 11, 3040-3047, 2011.
1809