Thin-Walled Structures 188 (2023) 110840

Contents lists available at ScienceDirect

Thin-Walled Structures
journal homepage: www.elsevier.com/locate/tws

Full length article

Free vibration analysis of a rotating skew plate by using the absolute nodal
coordinate formulation
Jing Zhang a , Xiaokang Du a , Yuanzhao Chen b , Xian Guo a , Liang Li a , Dingguo Zhang a ,∗
a
School of Science, Nanjing University of Science and Technology, Nanjing, 210094, Jiangsu, PR China
b
School of Mechanical and Automotive Engineering, Guangxi University of Science and Technology, Liuzhou, 545006, Guangxi, PR China

ARTICLE INFO ABSTRACT

Keywords: The skew angle is an important feature for the rotating plate, which has a significant influence on the modal
Absolute nodal coordinate formulation (ANCF) characteristics of the vibration. However, the variation of modal properties induced by the skew angle for a
Rotating skew plate rotating cantilever plate has seldom been studied, and most work has been focused on the static behaviors or
Free vibration
free vibration analysis of non-rotating skew plates in the existing literature. A novel dynamic model using a
Frequency loci veering
floating frame of reference (FFR) based on the absolute nodal coordinate formulation (ANCF) shell element
Critical speed
is proposed to investigate the vibration behaviors of a rotating skew plate. The convergence and validation
analyses are performed, and the results obtained by the present method are verified by comparing with those
in the available literature. A detailed parameters investigation on the effects of skew angle, aspect ratio, hub
radius and stagger angle on the variations of modal characteristics are demonstrated. The eigenvalue loci
veering phenomena along with the corresponding mode shape variations are presented and comprehensively
discussed. In addition, it is found that the skew plate may buckle due to the non-uniform centrifugal force,
which leads to the zero fundamental frequency by the critical rotational speed.

1. Introduction which are commonly encountered in the practical engineering field.

Therefore, attempts have been made to formulate the behaviors of
The rotating structures have a wide range of applications in various rotating blades with plate theory. Yoo et al. [7,8] derived a set of linear
engineering areas, such as blades of turbines, turbo engines and heli- equations for the modal analysis of a rotating plate. Fang et al. [9]
copters. The changes in dynamic behaviors and modal characteristics of studied the free vibration of rotating rectangle Mindlin plates with vari-
the flexible structure are significant when undergoing overall motions. able thickness via the Chebyshev–Ritz method. Li and Zhang et al. [10]
The results are attributed to the centrifugal inertia forces, which lead investigated the dynamic behaviors of a functionally graded rectangle
to the stretching of structures and the increase of their bending stiff- plate undergoing large overall motions. Most of the aforementioned
ness. The effect of Coriolis forces on the vibration couplings between
studies made the parameter investigation by considering the effect
different vibration modes is also considered in the research of rotating
of the aspect ratio, thickness ratio, hub radius, etc. Since the skew
structures. Therefore, analytical methods are desirable to identify the
angle is an important feature for the blade-like structure in practical
dynamic behaviors of these structures properly.
engineering, such as swept wings for aeroplanes, skew floor slabs and
In the early research work, beam models were extensively adopted
as rotating structures to analyze the dynamic response and vibration ship decks, the vibration characteristics of the skew plate were studied
characteristics. Southwell and Gough [1] studied the modal character- by using different numerical methods. The early studies [11–13] use
istics of rotating beams via the Southwell equation. Lee et al. [2,3] Rayleigh–Ritz method in the analysis of the vibration of the skew plate.
studied the free vibration of a beam rotating at a constant angular Durvasula [14] applied Galerkin’s method to investigate the frequencies
velocity. Yoo et al. [4] used the Rayleigh–Ritz method to investigate and mode shapes of clamped skew plates. Monforton [15] adopted the
the vibration of a rotating pre-twisted blade with a concentrated mass. finite element method to analyze the natural frequencies of clamped
Vallejo et al. [5,6] investigated the instability of a rotating beam rhombic plates. Garg et al. [16] studied the free vibration of isotropic,
by using the ANCF. They found that there exists a critical angular orthotropic and layered anisotropic composite laminated skew plates
velocity that made an unstable model of the beam. However, the beam through the iso-parametric finite element method. McGee and Kim [17]
models are not suitable for the accurate prediction of the high order proposed a modified Rayleigh–Ritz method, which applied B-spline
modal characteristics of the blade-type structure with low aspect ratio, functions into the procedure, to analyze the vibration and buckling

∗ Corresponding author.
E-mail address: zhangdg419@mail.njust.edu.cn (D. Zhang).

https://doi.org/10.1016/j.tws.2023.110840
Received 14 December 2022; Received in revised form 25 March 2023; Accepted 8 May 2023
Available online 23 May 2023
0263-8231/© 2023 Elsevier Ltd. All rights reserved.
J. Zhang, X. Du, Y. Chen et al. Thin-Walled Structures 188 (2023) 110840

of skew plates with both geometric and natural boundary conditions. vibration equation considering both centrifugal and Coriolis forces is
Kumar et al. [18–20] investigated the free vibration and buckling of established with the steady state of the rotating plate. The remainder of
laminated composite skew hyper shells using a 𝐶0 finite element for- the paper is organized as follows. The dynamic equations of a rotating
mulation baesd on high order shear deformation theory. As the reviews skew plate based on ANCF are derived in Section 2. In Section 3, the
demonstrate, most studies have focused on the vibration characteristics results from the presented method are compared with those from the
of the static skew plate with various boundary conditions, while the available literature to verify the accuracy of the method. In Section 4,
modal features with the rotating behavior are seldom investigated in the effects of the skew angle and the aspect ratio on the variations of
the existing literature. modal characteristics are performed in detail. In Section 5, the rotation-
To obtain the accurate dynamic behaviors of a flexible skew plate induced buckling is discussed and several parameters investigation of
subjected to overall rotation, the absolute nodal coordinate formulation the effects on variations of critical rotational speed is conducted. The
(ANCF) is utilized. The ANCF originally proposed by Shabana [21] results concluded from this research are presented in Section 6.
has been acknowledged as the non-incremental finite element method
(FEM) to deal with the dynamics of the flexible multibody systems un- 2. A skew plate model of ANCF
dergoing overall motions and large deformations. For the ANCF, global
position and the gradient vectors instead of infinitesimal rotations in In this study, two coordinate systems are established to describe the
conventional FEM, are used as the nodal coordinates to define the motion of the skew plate. One is the 𝑂−𝑋𝑌 𝑍 coordinate system, setting
element’s configuration. This leads to an exact description of rigid body as the inertial frame, and the other is the 𝑜 − 𝑥𝑦𝑧 coordinate system,
motion and the derived dynamic equation with constant mass matrix. which the origin lies at the root of the middle surface of the plate,
Various plate/shell elements based on ANCF have been proposed setting as the floating frame. As shown in Fig. 1, the position vector
for modeling plate and shell structures. Mikkola and Shabana [22] of an arbitrary point in the global inertial frame can be expressed as
proposed the four nodes fully parameterized ANCF plate element, 𝐫 = 𝐫𝑜 + 𝐫̃ (1)
which consists of the gradient vectors for all the coordinates parameters
and elicits a general description of element deformation including shear 𝐫̃ is the coordinate vector of the point 𝑃 in the rotating coordinate
deformation. Dimtrochenko and Pogorelov [23] proposed a gradient- frame of 𝑜 − 𝑥𝑦𝑧, and 𝐫𝑜 is the vector that indicates the original point
deficient thin plate element of ANCF. For the thin plate element, the of the inertial coordinate frame to the original point of the floating co-
configuration of the plate is characterized by the mid-surface only, in ordinate frame. The geometric parameters of the plate are represented
which a set of deficient gradient vectors is applied as nodal coordinates. by the span length 𝐿, the chord length 𝑊 and the thickness ℎ.
This thin plate element was verified with good accuracy and compu- Different configurations of a skew plate element of ANCF are shown
tational efficiency through the study by Dufva and Shabana [24]. Liu in Fig. 2. Based on the Kirchhoff–Love plate theory, the position of an
et al. [25] developed a thin cylindrical shell element of ANCF through arbitrary point on the surface 𝛤 𝑧 in the floating frame of the reference
the introduction of the angle between base vectors of a defined local configuration is defined as
frame. Yan et al. [26] adopted the approach to deal with the on-surface
𝐫̃ 0 (𝜉, 𝜂, 𝜁 ) = 𝐫̃ 𝑚0 + 𝜁 𝐧̃ 0 = 𝐒 (𝜉, 𝜂) 𝐞̃ 0 + 𝜁 𝐧̃ 0 (2)
slope discontinuity for modeling general curved structures based on the
gradient-deficient ANCF shell element. Valkeapaa et al. [27] proposed where the subscripts 𝑚 and 0 indicate that it is considered as on the
a shear deformable shell element based on the elastic middle surface mid-surface and in the initial configuration, respectively. 𝜉 and 𝜂 can be
approach using ANCF to study the method of locking remedies. Fur- used as the interpolation parameters of the mid-surface, 𝜁 denotes the
thermore, the strain split method [28] for the locking alleviation of distance between the mid-surface 𝛤 and 𝛤 𝑧 , and 𝐧̃ 0 is the unit normal
curved structures is applied to the ANCF shell element. vector of the mid-surface 𝛤 at point 𝑃 . 𝐞̃ 0 is the generalized coordinate
Recently, ANCF thin plate element with the non-local elasticity vector of the element with 36 degrees of freedom, which can be defined
theory is used to conduct a novel nonlinear model of single-layered as
graphene sheets subjected to large deformation [29]. Lan et al. [30] [ ( )T ( )T ( )T ( )T ]T
( )T 𝜕 𝐫̃ 0 𝜕 𝐫̃ 0 ( )T 𝜕 𝐫̃ 0 𝜕 𝐫̃ 0
developed a new ANCF thin plate tire model including the damp- 𝐞̃ 0 = 𝐫̃ 0 1 ... 𝐫̃ 0 4
ing property of the rubber to study the tire-ground contact process. 𝜕𝜉 1 𝜕𝜂 1 𝜕𝜉 4 𝜕𝜂 4
Yao et al. [31] studied the rotational vibration characteristics of the (3)
thin film antenna by using the ANCF with improved computational
efficiency. The transient dynamic behaviors and vibration analysis the subscripts 1, 2, 3, 4 denote the four nodes of the element. The
of a rotating plate with the stagger angle are investigated by Chen coordinate vector of node 1 can be written as
[ ( ) ( ) ]T
et al. [32]. The study indicates that the torsional deformation induced ( ) 𝜕 𝐫̃ 0 T 𝜕 𝐫̃ 0 T
by the stagger angle affects the modal characteristics of the rotating (̃𝐞0 )1 = 𝐫̃ 0 T1 (4)
𝜕𝜉 1 𝜕𝜂 1
plate significantly. Similar to the approach of Jiang et al. [33], the ( ) ( )
( ) 𝜕 𝐫̃ 0 𝜕 𝐫̃ 0
modal analysis is taken with the artificially centrifugal forces based on where 𝐫̃ 0 1 and , are the position vector and the slope
the ANCF. 𝜕𝜉 1 𝜕𝜂 1
vectors, respectively. 𝐒 is the shape function matrix of the element,
To the authors’ knowledge, there is a lack of research work devoted
which can be found in the work of Dufva and Shabana [21]. In
to the vibration analysis of a large overall rotating skew plate. In this
analogy, the position of the point 𝑃 in the current configuration can
paper, a novel structure model is developed to study the free vibration
be represented by substituting the initial generalized coordinate vector
characteristics of a rotating skew plate. Since the position vectors in the
𝐞̃ 0 in Eq. (2) with the current vector 𝐞̃ , which is
ANCF are described in the global coordinate system, which the large
overall rigid body motion and the flexible deformation displacement 𝐫̃ (𝜉, 𝜂, 𝜁 ) = 𝐒 (𝜉, 𝜂) 𝐞̃ + 𝜁 𝐧̃ (5)
cannot be separated directly, and to the dynamic model of rotational
According to the non-linear continuum mechanics, the Green strain
structures, the formulation of the centrifugal and the Coriolis forces
tensor at an arbitrary material point in the 𝑖th element can be defined
cannot be obtained directly by ANCF. Thus, the proposed model uses
as follows,
the nodal coordinates in the FFR to determine the displacement field
[( )T ] {[ }
and formulates the elastic energy based on the ANCF element, which is
𝑖 1 𝜕 𝐫̃ 𝑖 𝜕 𝐫̃ 𝑖 1 𝑖
( 𝑖 )−1 ]T [ 𝑖 ( 𝑖 )−1 ]
accurately described without the use of small deformation assumption. 𝐄 (𝜉, 𝜂, 𝜁 ) = − 𝐈 = 𝐉 𝐉0
𝐉 𝐉0
− 𝐈
2 𝜕𝐗𝑖 𝜕𝐗𝑖 2
The ANCF shell element is applied to derive the dynamic equation and
the steady state is obtained by solving the quasi-static equations. The (6)

2
J. Zhang, X. Du, Y. Chen et al. Thin-Walled Structures 188 (2023) 110840

Fig. 1. Configuration of a skew plate model.

Fig. 2. Different configurations of a skew plate element of ANCF.

where the 𝐗𝑖 represents the global position of the 𝑖th element in the where the unit vector along the 𝑍-axis is parallel to the unit normal
reference configuration. The covariant base tensor at the deformed and vector 𝐧̃ 0 , that is
reference configuration 𝐉𝑖 and 𝐉𝑖0 are defined as
𝑖
[ 𝑖 ] 𝐤 = 𝐧̃ 𝑖0 (12)
𝜕 𝐫̃ 𝑖 𝜕 𝐫̃ 𝜕 𝐫̃ 𝑖 𝜕 𝐫̃ 𝑖 [ ]
𝐉𝑖 = 𝑖 = = 𝐠𝑖1 𝐠𝑖2 𝐠𝑖3 (7)
𝜕𝐱 𝜕𝜁 𝜕𝜂 𝜕𝜉
[ 𝑖 ] and two other base unit vectors are defined by Gram–Schmidt orthog-
𝜕𝐗𝑖 𝜕 𝐫̃ 0 𝜕 𝐫̃ 0𝑖 𝜕 𝐫̃ 0𝑖 [ ]
𝐉𝑖0 = = = 𝐆𝑖1 𝐆𝑖2 𝐆𝑖3 (8) onalization as follows
𝜕𝐱 𝑖 𝜕𝜁 𝜕𝜂 𝜕𝜉
𝑖 𝐆𝑖1
Furthermore, the Green strain tensor in the Eq. (6) can be expressed as 𝐢 = (13)
|𝐆𝑖1 |
( )−T 𝑖 ( 𝑖 )−1
𝐄𝑖 (𝜉, 𝜂, 𝜁 ) = 𝐉𝑖0 ̃ 𝐉
𝐄 (9) 𝑖 𝑖 𝑖
0
𝐣 =𝐤 ×𝐢 (14)
where ̃𝑖
𝐄 is the Green–Lagrange strain tensor defined in the reference
(curvilinear) non-orthogonal coordinates, which can be written as The regularized Green strain tensor is then transformed with respect to
[( ) ( ) ]
̃ 𝑖 = 1 𝐉𝑖 T 𝐉𝑖 − 𝐉𝑖 T 𝐉𝑖
𝐄 (10)
an orthogonal coordinate frame in a matrix form as
0 0
2 [ 𝑖 ] [ 𝑖 ]
It is important to note that the strain tensor defined with the non- 𝜀11 𝜀𝑖12 ( 𝑖 )T 𝐸̃ 11 𝐸̃ 12
𝑖
= 𝐓 𝐓𝑖
orthonormal base vector cannot be directly correlated by the consti- 𝜀𝑖12 𝜀𝑖22 𝐸̃ 12
𝑖 𝐸̃ 22
𝑖

tutive equation. ([ 𝑖 𝑖 ] [ 𝑖 ])
1 𝑖T (𝐠1 )(𝐠1 ) (𝐠𝑖1 )(𝐠𝑖2 ) (𝐆1 )(𝐆𝑖1 ) (𝐆𝑖1 )(𝐆𝑖2 )
To this end, the Green–Lagrange strain tensor given by Eq. (10) = (𝐓 ) − 𝐓𝑖 (15)
should be transformed to which defined within the orthogonal coor- 2 (𝐠𝑖2 )(𝐠𝑖1 ) (𝐠𝑖2 )(𝐠𝑖2 ) (𝐆𝑖2 )(𝐆𝑖1 ) (𝐆𝑖2 )(𝐆𝑖2 )
dinate frame 𝐀𝑖 , and the three base vectors of the frame can be defined
where 𝐓𝑖 is the transformation matrix that can be written as
as
[ ] [( ) ]−1
𝑖 𝑖 𝑖 T
𝐀𝑖 = 𝐢 𝐣 𝐤 (11) 𝐓𝑖 = 𝐉𝑖0 𝐀𝑖 (16)

3
J. Zhang, X. Du, Y. Chen et al. Thin-Walled Structures 188 (2023) 110840

For the skew element with a skew angle 𝛼, the transformation matrix where 𝐃 is the coefficient matrix of Young modulus, which can be
𝐓𝑖 can be derived as expressed as

⎡1 tan 𝛼 ⎤
𝐸 ⎢
⎡ 1 𝜐 0 ⎤
𝑖 ⎢ 𝑙𝑖 𝑙𝑖 ⎥ 𝐃= 𝜐 1 0 ⎥ (28)
𝐓 =⎢ ⎥ (17) 1 − 𝜐2 ⎢ ⎥
⎢0 ⎣ 0 0 1−𝜐 ⎦
− 𝑖 ⎥⎦
sec 𝛼
⎣ 𝑤 𝐸, 𝜐 and ℎ are the Young modulus, the Poisson ratio and the thickness
where 𝑙𝑖 and 𝑤𝑖 are the edge length of arbitrary 𝑖th element. of the element, respectively.
For a thin plate element, the high order terms of 𝜁 coordinates in The elastic force vector can be calculated by
( )
the thickness direction can be neglected and the simplified Green strain 𝜕𝑈 T
𝐐𝑒 = = 𝐐𝑚
𝑒 + 𝐐𝑒
𝜅
(29)
tensor in Eq. (10) can be obtained as the form of 𝜕 𝐞̃
where the elastic forces associated with the mid-surface strain energy
̃𝑖 ≈ 𝐄
𝐄 ̃𝑖 +𝐄
̃𝑖 (18)
𝑚 𝑘 and bending strain energy 𝐐𝑚 𝜅
𝑒 , 𝐐𝑒 are derived as
̃𝑖 { ( 𝑚 )T ( 𝑚 )T
where the mid-surface strains 𝐄𝑚 are expressed as 𝐸ℎ 𝜕𝜀11 𝜕𝜀22
𝐐𝑚 = 𝜀 𝑚
+ 𝜀 𝑚
( )T ( )T 𝑒 ∫𝐴𝑒 1 − 𝜐2 11 𝜕 𝐞̃ 22 𝜕 𝐞̃
⎡ 𝜕 𝐫̃ 𝑖 𝜕 𝐫̃ 𝑚𝑖 𝑖
𝜕 𝐫̃ 𝑚0 𝑖 ⎤
𝜕 𝐫̃ 𝑚0
̃𝑖 = 1 ⎢ 𝑚 ⎥ [ ( 𝑚 )T ( 𝑚 )T ]
𝐄 − (19) 𝜕𝜀22 𝜕𝜀11
𝑚 2 ⎢ 𝜕𝐱𝑖 𝜕𝐱𝑖 𝜕𝐱𝑖 𝜕𝐱𝑖 ⎥ + 𝜐 𝜀𝑚 + 𝜀𝑚
⎣ ⎦ 11 22
(30)
𝜕𝐞̃ 𝜕 𝐞̃
̃ 𝑖 are expressed as
And the bending strains 𝐄 ( 𝑚 )T }
𝑘 𝜕𝜀12
( ) ( 𝑖 )T 𝑖 +4 (1 − 𝜐) 𝜀𝑚 12
d𝐴𝑒
⎡ 𝜕 𝐫̃ 𝑖 T 𝑖 𝜕 𝐧̃ 0 𝜕 𝐧̃ 0 ⎤ 𝜕 𝐞̃
̃𝑖 = 𝜁 ⎢
𝐄 𝑚 𝜕 𝐧̃
− ⎥ (20) { ( ) ( )
⎢ 𝜕𝐱𝑖 𝜕𝜅11 T 𝜕𝜅22 T
𝑘
𝜕𝐱 𝑖
𝜕𝐱 𝑖
𝜕𝐱𝑖 ⎥ 𝐐𝜅𝑒 = (
𝐸ℎ
) 𝜅11
⎣ ⎦ ∫𝐴𝑒 12 1 − 𝜐2
+ 𝜅22
𝜕 𝐞̃ 𝜕 𝐞̃
By using Weingarten’s formula in differential geometry, the compo- [ ( )T ( ) ]
𝜕𝜅22 𝜕𝜅11 T
nents of bending strain can be derived as + 𝜐 𝜅11 + 𝜅22 (31)
𝜕 𝐞̃ 𝜕 𝐞̃
( ) ( )
( 𝑖) 𝜕 𝐫̃ 𝑚𝑖 𝜕 𝐧̃ 𝑖 𝑖 𝜕𝐧
𝜕 𝐫̃ 𝑚0 ̃ 𝑖0 𝑖
𝜕 2 𝐫̃ 𝑚 𝜕 2 𝐫̃ 𝑚0
𝑖 ( ) }
̃ 𝜕𝜅12 T
𝐄 𝜅 𝑟𝑠 = 𝜁 − = −𝜁 ̃𝑖
𝐧 − 𝐧̃ 𝑖 +4 (1 − 𝜐) 𝜅12 d𝐴𝑒
𝜕𝜒 𝑟 𝜕𝜒 𝑠 𝜕𝜒 𝑟 𝜕𝜒 𝑠 𝜕𝜒 𝑟 𝜕𝜒 𝑠 𝜕𝜒 𝑟 𝜕𝜒 𝑠 0 𝜕 𝐞̃
(21) in which the derivatives of components of the strain vector 𝜺𝑚 and
curvature vector 𝜿 for the nodal coordinates 𝐞̃ are obtained as follows
where the 𝜒 𝑟 and 𝜒 𝑠 represent the parameters 𝜉 and 𝜂 with subscripts ( )
𝜕𝜀𝑚 𝑙𝑔 1
𝑟, 𝑠 = 1, 2. = (𝐓)𝑟𝑙 (𝐓)𝑠𝑔 𝐒T,𝜒 𝑟 𝐒,𝜒 𝑠 + 𝐒T,𝜒 𝑠 𝐒,𝜒 𝑟 𝐞̃ 𝑘 (32)
𝜕 𝑒̃𝑖 2 𝑘𝑖
The corresponding regularized strain components can be written in
matrix form as (
([ 𝑖 ] ⎛ ( ) ) ( ( ) ) ⎞
𝜕𝐍 𝜕𝐍
1 ( 𝑖 )T (𝐠1𝑚 )(𝐠𝑖1𝑚 ) (𝐠𝑖1𝑚 )(𝐠𝑖2𝑚 ) ⎜( 𝐫̃ 𝑓 ,𝜒 𝑟 𝜒 𝑠 𝜅𝑙𝑔 𝐍𝑓
𝑖
𝜺𝑚 = 𝐓 𝜕𝜅𝑙𝑔 ) 𝜕 𝐞̃ 𝑓 𝜕 𝐞̃ 𝑓 𝑖 ⎟
2 (𝐠𝑖2𝑚 )(𝐠𝑖1𝑚 ) (𝐠𝑖2𝑚 )(𝐠𝑖2𝑚 ) = (𝐓)𝑟𝑙 (𝐓)𝑠𝑔 ⎜ 𝐧𝑓 𝐒𝑓 ,𝜒 𝑟 𝜒 𝑠 𝑖 + 𝑖
− ⎟
𝜕 𝑒̃𝑖 ⎜ |𝐍| |𝐍|2 ⎟
[ 𝑖 ]) [ 𝑚 ] ⎜ ⎟
(𝐆1𝑚 )(𝐆1𝑚 ) (𝐆𝑖1𝑚 )(𝐆𝑖2𝑚 )
𝑖
𝑖
𝜀11 𝜀𝑚
12 ⎝ ⎠
− 𝐓 = 𝑚 (22)
(𝐆𝑖2𝑚 )(𝐆𝑖1𝑚 ) (𝐆𝑖2𝑚 )(𝐆𝑖2𝑚 ) 𝜀12 𝜀𝑚
22 (33)
The other subscripts in Eqs. (32) and (33) are 𝑙, 𝑔 = 1, 2, 𝑓 = 1, 2, 3
2 𝑖
2 𝑖
⎛⎡ 𝜕 𝐫̃ 𝑚 𝜕 2 𝐫̃ 𝑚𝑖⎤ ⎡ 𝜕 𝐫̃ 𝑚0 𝜕 2 𝐫̃ 𝑚0
𝑖
⎤⎞ and 𝑖, 𝑘 = 1, . . . , 36.
⎜⎢ 𝜕𝜉 2 ⋅ 𝐧̃ 𝑖 ⋅ 𝐧̃ 𝑖 ⎥ ⎢ 2
⋅ 𝐧̃ 𝑖0 ⋅ 𝐧̃ 𝑖0 ⎥⎟
( )T ⎜⎢ 𝜕𝜉𝜕𝜂 ⎥ ⎢ 𝜕𝜉 𝜕𝜉𝜕𝜂 ⎥⎟ 𝑖
𝜺𝑖𝜅 = −𝜁 𝐓𝑖 ⎜⎢ ⎥−⎢ ⎥⎟ 𝐓 3. Dynamic modeling of a rotational skew plate
⎜⎢ 𝜕 2 𝐫̃ 𝑖 𝜕 2 𝐫̃ 𝑚𝑖 ⎥ ⎢ 𝜕 2 𝐫̃ 𝑖 𝜕 2 𝐫̃ 𝑚0
𝑖 ⎥⎟
⎜⎢ 𝑚
⋅ 𝐧̃ 𝑖 ⋅ 𝐧̃ ⎥ ⎢
𝑖 𝑚0
⋅ 𝐧̃ 𝑖0 ⋅ 𝐧̃ 𝑖0 ⎥⎟
⎝⎣ 𝜕𝜉𝜕𝜂 𝜕𝜂 2 ⎦ ⎣ 𝜕𝜉𝜕𝜂 𝜕𝜂 2 ⎦⎠ A skew plate with a uniform mass distribution is molded as a
[ 0 0
] rotating skew blade which is clamped to a rigid disk with a radius 𝑅
𝜅11 − 𝜅11 𝜅12 − 𝜅12
= −𝜁
0 0
and mounted with a stagger angle 𝛹 , as shown in Fig. 3. The position
𝜅12 − 𝜅12 𝜅22 − 𝜅22
vector in the global inertial frame 𝑂 − 𝑋𝑌 𝑍 can be expressed as
(23) ( )
𝐫 = Φ𝜃𝜓 𝐫̃ ℎ𝑢𝑏 + 𝐫̃ (34)
For fitting the calculation of the linear constitutive relation, the 𝐫̃ ℎ𝑢𝑏 is the vector that indicates the rotation axis of the rigid hub to the
regularized strain tensor can be expressed as the vector of Voigt form original point of the rotating coordinate frame 𝑜 − 𝑥𝑦𝑧, which can be
[ ]T expressed as
𝜺 = 𝜀11 𝜀22 2𝜀12 = 𝜺𝑚 − 𝜁𝜿 (24)
[ ]T
𝐫̃ ℎ𝑢𝑏 = 𝑅 0 0 (35)
in which
[ ]T where 𝑅 is the radius of the hub rigid. The coordinate transformation
𝜺 𝑚 = 𝜀𝑚11
𝜀𝑚
22
2𝜀𝑚
12
(25) matrix Φ𝜃𝜓 is written as

[ 0 0 0 )
]T Φ𝜃𝜓 = Φ𝑌 ,𝜃 Φ𝑋,𝜓 (36)
𝜿= 𝜅11 − 𝜅11 𝜅22 − 𝜅22 2(𝜅12 − 𝜅12 (26)
in which the components matrixes are defined as the following expres-
Based on the Kirchhoff–Love plate theory, the strain energy of the sion
element can be written as
( ) ⎡ cos 𝜃 0 sin 𝜃 ⎤
1 ℎ ℎ2 T Φ𝑌 ,𝜃 = ⎢ 0 1 0 ⎥ (37)
𝑈= 𝜺T 𝐃𝜺d𝑉𝑒 = (𝜺𝑚 )T 𝐃 (𝜺𝑚 ) + 𝜿 𝐃𝜿 d𝐴𝑒 (27) ⎢ ⎥
2 ∫𝑉𝑒 2 ∫ 𝐴𝑒 12 ⎣− sin 𝜃 0 cos 𝜃 ⎦

4
J. Zhang, X. Du, Y. Chen et al. Thin-Walled Structures 188 (2023) 110840

∑𝑚
𝐅𝑐 = − 𝐒T ΦT𝜃𝜓 Ω
̂ ΩΦ
̂ 𝜃𝜓 𝐫 𝑖 d𝑉𝑒 (47)
∫
𝑖=1 𝑉𝑒

𝐪 is the generalized coordinate vector, which can be expressed as

[ ( ) ( 𝑖 )T ]T
T
𝐪= 𝐞̃ 1 ⋯ 𝐞̃ ⋯ (̃𝐞𝑚 )T (48)

Φ(𝐪, 𝑡) is the system constraint equations, Φ𝐪 is the derivative matrix of

constraint equations with respect to the 𝐪, and 𝝀 is the Lagrange multi-
plier. The term underline represents the effect of angular acceleration
to the system, which can be neglected when the plate rotates with a
steady angular velocity. In the dynamic model of a rotating skew plate,
as the rotation axis paralleling to the 𝑌 -axis of the global inertial frame
𝑂 − 𝑋𝑌 𝑍, the angular velocity vector Ω can be written as
[ ]T
Ω = 0 𝜃̇ 0 (49)
̂ in Eq. (47) is
And the corresponding matrix Ω
Fig. 3. Configuration of a rotational skew plate model. 𝜃̇
⎡ 0 0 ⎤
̂ =⎢ 0
Ω 0 0 ⎥ (50)
⎢ ⎥
⎣ −𝜃̇ 0 0 ⎦
⎡1 0 0 ⎤ By substituting Eqs. (34) and (50) into Eq. (47), the centrifugal
Φ𝑋,𝜓 = ⎢0 cos 𝜓 − sin 𝜓 ⎥ (38) forces vector is re-expressed as follows,
⎢ ⎥
⎣0 sin 𝜓 cos 𝜓 ⎦ ∑
𝑚
( )
𝐅𝑐 = −𝜃̇ 2 𝐒T 𝐒̃𝐞𝑖 + 𝐫̃ ℎ𝑢𝑏 d𝑉𝑒 (51)
The kinetic energy of each element can be obtained as ∫
𝑖=1 𝑉𝑒
1 𝐒 is a 3 × 36 matrix, in which the components of the first and third
𝑇 = 𝜌 𝐫̇ T 𝐫d𝑉
̇ 𝑒 (39)
2 ∫𝑉𝑒 rows are equal to the corresponding components of the shape function
in which the velocity vector 𝐫̇ in the global inertial frame can be derived 𝐒, while the other components are equal to zero, that is
by differentiation of the position vector 𝐫 in Eq. (34) with respect to [ ]T
𝐒 = 𝐒1 T 0 𝐒3 T (52)
time, that is
Since the constraint equations Φ(𝐪, 𝑡) are linear in the rotating
̂ + 𝛷𝜃𝜓 𝐒𝐞̃̇
𝐫̇ = 𝜴 × 𝐫 + Φ𝜃𝜓 𝐫̃̇ = 𝛺𝑟 (40) frame, it is possible and easy to get rid of the constraints by using the
Ω is the angular velocity of the rotating hub rigid defined in the transformation
𝑂 − 𝑋𝑌 𝑍 coordinate and Ω̂ is a 3 by 3 matrix defined as
𝐪 = 𝐁𝐪∗ (53)
⎡ 0 −𝛺3 𝛺2 ⎤ where 𝐪∗
is the independent generalized coordinate vector. When
̂ = ⎢ 𝛺3
Ω 0 −𝛺1 ⎥ (41)
⎢ ⎥ considering the situation of rotating plate with steady angular velocity,
⎣ −𝛺2 𝛺1 0 ⎦ the Eq. (43) can be simplified as
where 𝛺1 , 𝛺2 and 𝛺3 are the three components of the angular velocity ∗
𝐌∗ 𝐪̈ ∗ + 𝐂∗ 𝐪̇ ∗ + 𝐐int = 𝐅𝑐 ∗ (54)
vector Ω.
By substituting the Eq. (40) into Eq. (39), the expression of the where
kinetic energy can be expanded as follows, ∗
𝐌∗ = 𝐁T 𝐌𝐁, 𝐂∗ = 𝐁T 𝐂𝐁, 𝐐int = 𝐁T 𝐐int , 𝐅𝑐 ∗ = 𝐁T 𝐅𝑐 (55)
[ ( ) ]
1 ̂ T Φ𝜃𝜓 𝐒𝐞̃̇ + 𝐞̃̇ T 𝐒T 𝐒 𝐞̃̇ d𝑉𝑒
𝑇 = 𝜌 𝐫 T ΦT𝜃𝜓 Ω
̂ T ΩΦ
̂ 𝜃𝜓 𝐫 + 2𝐫 T ΦT Ω
𝜃𝜓
In order to obtain the eigenvalue equations for vibration analysis of
2 ∫𝑉𝑒 a rotating skewed plate, the vibration eigenvalue equation should be
(42) performed linearized with the steady state of the generalized coordi-
nate vector, in which the effect of non-zero centrifugal force is taken
Based on the ANCF, the final dynamic equations for a system of a into consideration. By neglecting the angular acceleration term, the
rotating skew plate can be expressed as differential–algebraic equations perturbation form of the dynamic equation can be written as
as follows
𝐌∗ 𝛿 𝐪̈ ∗ + 𝐂∗ 𝛿 𝐪̇ ∗ + 𝐊𝑇 ∗ 𝛿𝐪∗ = 𝟎 (56)
⎧ ∑𝑚
⎪ 𝐌𝐪̈ + 𝐂𝐪̇ + ΦT𝐪 𝝀 + 𝐐int + ̂̇ 𝜃𝜓 𝐫 𝑖 = 𝐅𝑐
𝐒T ΦT𝜃𝜓 ΩΦ ∗
⎨ (43) 𝐊𝑇 is the tangential stiffness matrix which can be defined as
⎪
𝑖=1 ( )
Φ (𝐪, 𝑡) = 𝟎 ∗
⎩ 𝜕 𝐐int − 𝐅𝑐 ∗
𝐊𝑇 ∗ = (57)
where the constant mass matrix 𝐌, the Coriolis forces matrix 𝐂, the 𝜕𝐪∗
generalized elastic forces vector 𝐐𝑖𝑛𝑡 and the centrifugal forces vector To obtain the generalized coordinate vector of static equilibrium
𝐅𝑐 can be defined as follows, configuration, the static equilibrium equation is solved as follows
∑
𝑚
∗
𝐌= 𝜌 𝐒T 𝐒d𝑉𝑒 (44) 𝐐int − 𝐅𝑐 ∗ = 𝟎 (58)
∫𝑉𝑒
𝑖=1
And Eq. (58) can be solved by using Newton–Raphson method. Since
∑
𝑚
the dynamic equation in Eq. (56) cannot be decoupled with the physical
T
𝐂= 2𝜌 𝐒 ΦT𝜃𝜓 ΩΦ
̂ 𝜃𝜓 𝐒d𝑉𝑒 (45)
∫𝑉𝑒 and the real mode coordinates, the description of state space is applied.
𝑖=1
∑
𝑚 The state variable 𝐮 is defined as
𝐐int = 𝐐𝑒 (46) [ ]T
𝐮 (𝑡) = (𝐪∗ )T (𝐪̇ ∗ )T (59)
𝑖=1

5
J. Zhang, X. Du, Y. Chen et al. Thin-Walled Structures 188 (2023) 110840

Table 1
Comparison of lowest five natural frequencies.
Dimensionless Vibration Yoo’s ANSYS Present method errors (%)
rotational speed model method [34]
2 × 2 4 × 4 8 × 8 12 × 12 16 × 16
1 3.5963 3.5798 3.5697 3.5779 3.5796 3.5799 3.58 0.0028
2 8.5507 8.4987 8.5754 8.5366 8.525 8.5226 8.5218 0.0094
𝛾=2 3 21.865 21.724 22.2723 21.9975 21.8106 21.7712 21.7571 0.065
4 27.384 27.193 26.3483 27.0387 27.1909 27.2202 27.2304 0.037
5 31.477 31.217 30.6053 31.2742 31.3172 31.3134 31.3107 0.0086
1 5.0491 5.0712 5.3338 5.0731 5.0312 5.0282 5.0279 0.006
2 9.0322 8.8365 9.1286 8.8738 8.8091 8.8 8.7973 0.031
𝛾 = 10 3 26.761 26.359 26.3921 26.3543 26.2962 26.2884 26.2861 0.0087
4 32.35 32.581 32.5384 32.5064 32.4066 32.3915 32.3872 0.013
5 39.078 38.933 38.4899 38.8356 38.7376 38.716 38.7087 0.019

By substituting Eq. (59) into Eq. (56), one can obtain the first order
linear differential equation of state space representation
[ ∗ ] [ ]
𝐂 𝐌∗ 𝐊𝑇 ∗ 𝟎
∗ 𝛿 𝐮̇ + ∗ 𝛿𝐮 = 𝟎 (60)
𝐌 𝟎 𝟎 −𝐌
The solution of Eq. (60) can be written as

𝛿𝐮 = 𝑒j𝜔𝑡 𝐙 (61)

𝐙 is a constant complex vector that identifies the deflection of

the mode shape. Substituting Eq. (61) into Eq. (60), the eigenvalue
equation is expressed as
[ ∗ ] [ ]
𝐂 𝐌∗ 𝐊𝑇 ∗ 𝟎
𝜔j 𝐙 = − 𝐙 (62)
𝐌∗ 𝟎 𝟎 −𝐌∗
The natural frequencies and associated mode shapes of the plate can
be obtained by solving the above equation. All the numerical results are
presented as dimensionless parameters, which are defined as follows
√
𝐿 ℎ 𝑅 ∗ ∗ 𝜌ℎ𝐿4
𝛿= , 𝜂 = , 𝜎 = , 𝛾 = 𝛺𝑇 , 𝑇 = , 𝜔̄ = 𝜔𝑇 ∗ Fig. 4. Transient deformations of rotating plate with different skew angles.
𝑊 𝐿 𝐿 𝐷
(63)

where
𝐸ℎ3
𝐷= ( ) (64)
12 1 − 𝜐2
𝛿, 𝜂, 𝜎 and 𝛹 denote the aspect ratio, thickness ratio, hub radius ratio
and the stagger angle of the skew plate, respectively. 𝛾 and 𝜔 represent
the dimensionless rotational speed and natural frequency.

4. Convergence and verification study

The convergence study of the present modeling method is verified

by comparing the first five dimensionless frequencies obtained by the
developed model with those in Yoo and Pierre [34] and the results from
ANSYS (see Table 1).
The dimensionless parameters are adopted as 𝛿 = 1, 𝜂 = 0.01,
𝜎 = 0. It is found that when the number of elements in this numerical
calculation is more than 12 × 12, results with enough accuracy have
been observed, which the maximal difference between 12 × 12 and
16 × 16 meshes is less than 0.1%. Hence, a 12 × 12 mesh of elements
is employed in the following study. Fig. 5. First six dimensionless frequencies of the plate versus dimensionless rotational
To prove the reliability of the present method for the plate with speed with 𝛼 = 0◦ and 𝛼 = 15◦ .
skew angle, a comparison of the first five dimensionless frequencies of
cantilever skew plate with the results in the available literature is given
in Table 2. The dimensionless parameters and material properties are
the dimensionless frequencies obtained from the present method are
selected as 𝛿 = 1, 𝜂 = 0.01, 𝜎 = 0, 𝛹 = 0◦ , 𝜌=2750 kg∕m3 , 𝐸 = 7×1010
slightly higher compared with those results of previously literature.
N∕m2 and 𝜐 = 0.3.
This could be addressed to the fact that the influence of the corner stress
It can be observed that for the skew angle 𝛼 = 15◦ , 30◦ , 45◦ , there is
singularities becomes significant with the large skew angle.
a good agreement between the results from the present method and
those from the Ritz method [35,37,38] and the conventional finite To predict transient dynamics of the rotating plates, a numerical
element [36] solutions. Moreover, when the skew angle 𝛼 = 60◦ , case with known angular speed is considered. The geometric and

6
J. Zhang, X. Du, Y. Chen et al. Thin-Walled Structures 188 (2023) 110840

Table 2
Comparison of lowest five natural frequencies.
Dimensionless Vibration K. M. Liew’s W. H. Liu’s L. Zhou’s Oliver G.’s Present
rotational speed model method [35] method [36] method [37] method method
[38]
1 3.5837 3.5884 —— 3.5813 3.5834
2 8.7918 8.713 —— 8.6968 8.6999
𝛼 = 15◦ 3 22.2372 22.2818 —— 22.229 22.277
4 27.6053 —— —— 26.333 26.3024
5 34.1824 —— —— 33.961 33.8786

1 3.9311 3.9383 3.9268 3.9278 3.9315

2 9.4126 9.4618 9.4096 9.4098 9.4177
◦
𝛼 = 30 3 25.2968 25.4436 25.28 25.286 25.4091
4 25.9344 —— 25.928 25.931 25.9113
5 41.3408 —— 41.332 41.331 41.3713

1 4.5114 4.5469 —— 4.5051 4.5306

2 11.2553 11.4237 —— 11.247 11.2498
𝛼 = 45◦ 3 27.0003 27.4364 —— 26.967 27.0692
4 31.565 —— —— 31.503 31.9204
5 50.715 —— —— 50. 710 50.485

1 5.2763 —— 5.2515 5.2427 5.2959

2 16.0973 —— 16.017 16.022 16.0097
𝛼 = 60◦ 3 30.658 —— 30.408 30.357 30.7281
4 45.6035 —— 45.348 45.29 45.8033
5 59.3242 —— 59.094 59.035 59.3037

Fig. 6. Variation of nodal lines of the third mode shape with an increase of rotational speed.

material parameters are given as 𝐿 = 1.0 m, 𝑊 = 0.5 m, ℎ = 0.0025 m, 5. Free vibration of the rotating skew plate
𝜌 = 2750 kg∕m3 , 𝐸 = 7 × 1010 N∕m2 and 𝜐 = 0.3. The angular speed of
the plate is given as
In this section, the effects of skew angle and aspect ratio on the
( ) free vibration of the rotating skew plate are studied. In the following
⎧ 𝛺𝑠 𝛺𝑠 2𝜋
⎪ 𝑡− sin 𝑡 , 0 ≤ 𝑡 ≤ 𝑇0 discussion, if not declared specially, the parameters of the skew plate
𝛺 = ⎨ 𝑇𝑠 2𝜋 𝑇𝑠 (65)
⎪𝛺 , are set as follows: 𝛿 = 1, 𝜂 = 0.01, 𝜎 = 0 and 𝛹 = 0◦ .
⎩ 𝑠 𝑡 > 𝑇0
Fig. 5 shows the corresponding first six dimensionless frequencies
of the rotating skew plate versus the dimensionless rotational speed.
Fig. 4 shows the variation of the deformations 𝑢𝑧 of the corner points It is observed that all the frequencies increase as the angular speed
on the free edges of the plate with time. The steady state speed 𝛺𝑠 and increase for the case of skew angle 𝛼 = 0◦ , while the fundamental
the acceleration time 𝑇𝑠 are set as 50 rad∕s and 5 s. It is obviously seen and third frequencies firstly increase and then turn slightly decrease
from the figure that the rotating plate with non-zero skew angles will for 𝛼 = 15◦ . Moreover, the second frequency of the rotating skew plate
have a large lateral deformation when the plate rotates at the steady keeps decreasing as the rotational speed increases. Similar to the case
state. It can be speculated that the skew angles on transient dynamics of the rotating rectangle plate, the frequency loci veering phenomenon
of rotating plates will influence the steady configurations significantly. can be observed obviously for the plate with a non-zero skew angle,

7
J. Zhang, X. Du, Y. Chen et al. Thin-Walled Structures 188 (2023) 110840

Fig. 7. Variation of nodal lines of the fourth mode shape with an increase of rotational speed.

Fig. 8. Variation of nodal lines of the third and fourth mode shapes versus dimensional
rotational speed: (a) the third mode shape; (b) the fourth mode shape.

where the third and fourth frequencies approach each other around the
velocity of 𝛾 = 7 and veer away with the increase of 𝛾.
Figs. 6–7 display the nodal line patterns of the third and fourth
frequency loci veering with various dimensionless rotational speeds,
respectively. As shown in Figs. 6 and 7, when the skew plate rotates
Fig. 9. First six dimensionless frequencies versus skew angle for the plate with 𝛾 = 2.
at a velocity of 𝛾 = 6, the third and fourth modes are the second
span-wise bending dominant and the coupled dominant of chordwise
bending mode and torsional mode. However, when the dimensionless
rotational speed is 𝛾 = 8, the third mode becomes the chord-wise fourth mode when the skew angle is around 𝛼 = 30◦ . Fig. 10 shows the
bending dominant and the fourth mode becomes the second span-wise corresponding mode shapes in detail.
bending dominant. Fig. 8 shows the variation of the nodal line patterns Fig. 10(a–c) show the third mode shapes for skew angle 𝛼 = 15◦ , 30◦
from Figs. 6 and 7. It can be seen that the mode loci veering occurs and 45◦ , and Fig. 10(d–f) show the fourth mode shapes for the same
continuously, in which all the nodal lines of the two mode shapes pass variation of skew angle when 𝛾 = 2. It can be observed that the third
through two fixed points in the plate. mode shape for represents the second spanwise bending mode, and the
Since the skew angle has an influence on the frequencies of the fourth mode shape represents the chordwise bending mode. While the
rotating plate, further study on the effect of skew angle is conducted third mode shape for 𝛼 = 15◦ represents the chordwise bending mode,
herein. Fig. 9 indicates the first six dimensionless frequencies of the and the fourth mode shape represents the second spanwise bending
plate versus the variation of skew angle 𝛼 with dimensionless rotational mode. The mode exchange between the second spanwise bending mode
speed 𝛾 = 2. and chordwise bending mode occurs with the increase of skew angle.
It is shown that the frequency veering phenomenon happens be- The variation of the nodal lines for the first six frequencies versus
tween the third and fourth modes. The third mode exchanges with the the skew angle is displayed in Figs. 11–14.

8
J. Zhang, X. Du, Y. Chen et al. Thin-Walled Structures 188 (2023) 110840

Fig. 10. Mode shapes of the skew plate with 𝛾 = 2.

Fig. 11. First six nodal line patterns for the skew plate with 𝛼 = 15◦ , 𝛾 = 2.

It can be found that when the plate has a large skew angle, the torsional mode develops the nodal line shifting from the center to the
fundamental mode shows the coupling between the spanwise bending bottom of the plate with the increase of the skew angle. Moreover, the
and torsional vibration, where the nodal line occurs in the first mode two nodal lines of the fifth mode, which is dominated by the second
nodal pattern of the plate. Furthermore, the nodal pattern of the first torsional mode, separate further from each other with the increase of

9
J. Zhang, X. Du, Y. Chen et al. Thin-Walled Structures 188 (2023) 110840

Fig. 12. First six nodal line patterns for the skew plate with 𝛼 = 30◦ , 𝛾 = 2.

Fig. 13. First six nodal line patterns for the skew plate with 𝛼 = 45◦ , 𝛾 = 2.

the skew angle. It indicates that there exists a strong coupling between turn slightly decrease as the aspect ratio increase then exchanges mode
chordwise bending and the second spanwise bending mode of the with the second frequency at 𝛿 = 3. Analogously, the mode exchanges
rotating plate with a large skew angle. between fifth and sixth frequencies with 𝛾 = 10 occurs at the aspect
The effect of the aspect ratio 𝛿 on the variation of the frequencies ratio of around 𝛿 = 3.5. The corresponding mode shapes are shown
is shown in Fig. 15. It can be found that the variations of the first in Figs. 16–17. It can be observed that the second and third modes
two frequencies with rotation speed 𝛾 = 2 and 𝛾 = 10 are similar for 𝛿 = 2.5 with 𝛾 = 2 are the first torsional-dominant and spanwise
to the increase of the aspect ratio. As the aspect ratio increases, the bending-dominant modes. However, when the aspect ratio is 𝛿 = 3.5,
fundamental frequency decreases while the second and fourth frequen- the second and third modes shift to each other compared to the 𝛿 = 2.5.
cies increase. The third frequency with 𝛾 = 2 firstly increases and A similar mode veering phenomenon occurs between the aspect ratio

10
J. Zhang, X. Du, Y. Chen et al. Thin-Walled Structures 188 (2023) 110840

Fig. 14. First six nodal line patterns for the skew plate with 𝛼 = 60◦ , 𝛾 = 2.

Fig. 15. First six dimensionless frequencies versus aspect ratio for the plate with skew angle 𝛼 = 15◦ .

𝛿 = 3 and 𝛿 = 4 when 𝛾 = 10, where the fifth and sixth modes are 6. Critical rotational speed for fundamental frequencies
dominated by the second torsional and the first torsional modes for
𝛿 = 3 and exchange to each other for 𝛿 = 4. As mentioned in Section 4, the fundamental frequency of the rotat-
Figs. 18–21 display the first six nodal line patterns for the skew plate ing skew plate will decrease with a high rotational speed, which can
with 𝛿 = 0.5, 1, 2.5, 4, respectively. It can be observed that the node line be seen in Fig. 5. Therefore, studies for the fundamental frequency are
patterns for the fundamental mode show the coupling between the first taken with several parameters investigation.
spanwise bending and torsional vibration. The first two mode shapes It is seen from Fig. 22 that the fundamental dimensionless fre-
are the same with the variation of the increased aspect ratio, while the quencies of the rotating skew plates decrease with increasing rotation
third to the sixth mode shapes are complicated and different with these speed 𝛾 when skew angle 𝛼 is non-zero. At a certain rotation speed,
aspect ratios. Moreover, as the aspect ratio increases, the mode shape of the frequencies reduce to zeros, which indicates that the buckling of
chordwise bending vibration shifts from a lower frequency to a higher the plate occurs. This is mostly attributed to the non-uniform cen-
one, which can be attributed to the plate’s effective length in spanwise trifugal forces of the rotating skew plate, which induce the coupled
bending increase and the relative extensional stiffness decrease, causing longitudinal extensional and lateral–torsional deflection of the plate.
the less low modes to have the mode patterns with longitudinal nodal The equilibrium configuration is instability when the dimensionless
lines. rotational speed increase to the critical speed and the determinant of

11
J. Zhang, X. Du, Y. Chen et al. Thin-Walled Structures 188 (2023) 110840

Fig. 16. Second and third nodal line patterns for the skew plate with 𝛿 = 2.5, 3, 3.5, 𝛼 = 15◦ , 𝛾 = 2.

Fig. 17. Fifth and sixth nodal line patterns for the skew plate with 𝛿 = 3, 3.5, 4, 𝛼 = 15◦ , 𝛾 = 10.

the stiffness matrix 𝐊𝑇 ∗ in Eq. (57) is equal to zero, which causes the length of the plate in spanwise increase with the increasing of aspect
zero fundamental frequency. It indicates that critical rotation speed ratio when 𝛿 > 1, which leads to the increasing of buckling load induced
decreases with the increase of skew angle 𝛼. For the case of non- by centrifugal forces.
rotating skew plate, the effective length in spanwise decreases with Variation of fundamental frequencies versus rotational speed for a
the increase of skew angle, which lead to the increase of the relative rotating skew plate with various hub radius ratios has been given in
bending stiffness. As the increase of rotational speed, the equilibrium Fig. 24. Buckling is observed when 𝛼 = 30◦ for all of the hub radius ra-
configuration of the plate is strong influenced by the non-uniform tios considered. It is observed that all the fundamental frequencies first
centrifugal forces and the instability state is more likely to occur. increase with the increasing rotational speed and then turn decrease.
The effect of the aspect ratio 𝛿 on the critical rotational speed is In addition, the peaks of the frequencies increase with the increase of
presented in Fig. 23. It can be seen that when the aspect ratio 𝛿 > 1, the the hub radius ratios. This is because the centrifugal forces increase
critical speed decreases with the increase of the aspect ratio. However, for the rotational speed increases as well as the hub radius. It should
when the aspect ratio 𝛿 = 0.5, the critical speed is lower than that of be noted that the critical rotational speed increase as the hub radius
the plate with 𝛿 = 1. The reason can be attributed to that the effective ratio increase when 𝜎 < 1, and on the contrary, the critical rotational

12
J. Zhang, X. Du, Y. Chen et al. Thin-Walled Structures 188 (2023) 110840

Fig. 18. Fifth and sixth nodal line patterns for the skew plate with 𝛼 = 15◦ , 𝛾 = 10, 𝛿 = 0.5.

Fig. 19. First six nodal line patterns for the skew plate with 𝛼 = 15◦ , 𝛾 = 10, 𝛿 = 1.

Fig. 20. First six nodal line patterns for the skew plate with 𝛼 = 15◦ , 𝛾 = 10, 𝛿 = 2.5.

speed turns to decrease with the increase of the hub radius ratio when by centrifugal forces is more significant than the tension effect of the
𝜎 > 1. It can be addressed to that when 𝜎 > 1 the buckling load induced rotation with the high rotational speed.

13
J. Zhang, X. Du, Y. Chen et al. Thin-Walled Structures 188 (2023) 110840

Fig. 21. First six nodal line patterns for the skew plate with 𝛼 = 15◦ , 𝛾 = 10, 𝛿 = 4.

Fig. 22. The fundamental frequency of the rotating skew plate versus rotational speed
with various skew angle 𝛼.
Fig. 24. The fundamental frequency of the rotating skew plate versus rotational speed
with various hub radius ratio 𝜎 for 𝛼 = 30◦ .

Fig. 25 shows the variation of the fundamental frequencies of the

rotating skew plate with different stagger angles versus the rotational
speed. As shown in the figure, all the fundamental frequencies of
the plate with non-zero stagger angle increase as the rotational speed
increases, where the phenomena of buckling are eliminated. And the
rate of increase becomes higher as the stagger angle 𝛹 increases. That is
because the centrifugal forces have a large lateral component with the
non-zero stagger angle, which leads to the stable torsional-dominant
steady state of the skew plate with the increasing rotational speed. The
fundamental frequencies always increase as the stagger angle increases
up to 60◦ , the main reason is that with the torsional-dominant steady
state, the rotational tension effect significantly when the stagger angle
increases.

7. Conclusion

In this paper, a novel dynamic model was developed to study the

Fig. 23. The fundamental frequency of the rotating skew plate versus rotational speed
with various aspect ratio 𝛿 for 𝛼 = 30◦ .
vibration characteristics of a rotating skew plate. The displacement
field is determined using the FFR and the elastic energy is formu-
lated accurately based on the ANCF shell element theory. Both of the
centrifugal and the Coriolis forces were considered in the established
dynamic equations. The accuracy of the present method and calculation

14
J. Zhang, X. Du, Y. Chen et al. Thin-Walled Structures 188 (2023) 110840

Acknowledgments

This research is funded by the grants from the National Natu-

ral Science Foundation of China (Project Nos. 12232012, 12072159,
12102191, and 12202110) and the Fundamental Research Funds for
the Central Universities, China (Project No. 30922010314).

References

[1] R. Southwell, F. Gough, The Free Transverse Vibration of Airscrew Blades, British
ARC Reports and Memoranda, No. 766, 1921.
[2] S.Y. Lee, S.M. Lin, C.T. Wu, Free vibration of a rotating non-uniform beam with
arbitrary pretwist, an elastically restrained root and a tip mass, J. Sound Vib.
273 (3) (2004) 477–492.
[3] S.M. Lin, The instability and vibration of rotating beams with arbitrary pretwist
and an elastically restrained root, J. Appl. Mech. 68 (6) (2001) 844–853.
[4] H.H. Yoo, J.Y. Kwak, J. Chung, Vibration analysis of rotating pre-twisted blades
with a concentrated mass, J. Sound Vib. 240 (5) (2001) 891–908.
[5] D. Garcia-Vallejo, H. Sugiyama, A.A. Shabana, Finite element analysis of the
geometric stiffening effect Part 1: a correction in the floating frame of reference
formulation, Proc. Inst. Mech. Eng. K 219 (2) (2005) 187–202.
[6] D. García-Vallejo, H. Sugiyama, A.A. Shabana, Finite element analysis of the
Fig. 25. The fundamental frequency of the rotating skew plate versus rotational speed geometric stiffening effect. Part 2: non-linear elasticity, Proc. Inst. Mech. Eng. K
with various stagger angle 𝛹 for 𝛼 = 30◦ . 219 (2) (2005) 203–211.
[7] H.H. Yoo, S.K. Kim, D.J. Inman, Modal analysis of rotating composite cantilever
plates, J. Sound Vib. 258 (2) (2002) 233–246.
[8] H.H. Yoo, S.K. Kim, Free vibration analysis of rotating cantilever plates, AIAA
results is confirmed by comparing with those in the available literature. J. 40 (11) (2002) 2188–2196.
The following conclusions are concluded: [9] J.S. Fang, D. Zhou, Free vibration analysis of rotating Mindlin plates with
variable thickness, Int. J. Struct. Stab. Dyn. 17 (04) (2017) 1750046.
1. The effects of the rotational speed on the frequencies and the [10] L. Li, D.G. Zhang, Free vibration analysis of rotating functionally graded
corresponding mode shapes of the skew plate were discussed. rectangular plates, Compos. Struct. 136 (2016) 493–504.
[11] M.V. Barton, Vibration of rectangular and skew cantilever plates, J. Appl. Mech.
The frequency loci veering can be observed and the nodal lines 18 (2) (1951) 129–134.
patterns of mode shape switch to each other continuously. [12] R.K. Kaul, V. Cadambe, The natural frequencies of thin skew plates, Aeronaut.
2. The effects of the skew angle and the aspect ratio on the varia- Q. 7 (4) (1956) 337–352.
tions of modal characteristics were investigated. The frequency [13] M. Hasegawa, Vibration of clamped parallelogrammic isotropic flat plates, J.
Aeronaut. Sci. 24 (2) (1957) 145–146.
loci veering with associated modes switching can be observed as
[14] S. Durvasula, Natural frequencies and modes of clamped skew plates, AIAA J. 7
the variation of frequencies with varying skew angles and aspect (6) (1969) 1164–1167.
ratios. The corresponding nodal lines patterns of mode shapes [15] G.R. Monforton, Finite Element Displacement and Vibration Analysis of Skew
indicated that the increase of skew angle will lead to the strong Plates, Report18, 1968, Division of Solid Mechanics.
coupling between the spanwise bending and torsional modes, as [16] A.K. Garg, R.K. Khare, T. Kant, Free vibration of skew fiber-reinforced composite
and sandwich laminates using a shear deformable finite element model, J. Sandw.
well as the coupling between the chordwise bending and second
Struct. Mater. 8 (1) (2006) 33–53.
spanwise bending modes. [17] O.G. McGee III, Flexural vibrations of clamped-free rhombic plates with corner
3. The investigation of the effects of the skew angle, aspect ra- stress singularities. Part I: review of research, J. Vib. Control 21 (13) (2015)
tio, hub radius and stagger angle on the critical speed was 2583–2603.
conducted. The results demonstrated that the critical speed de- [18] A. Kumar, C. Anupam, K. Mrunal, Analysis of laminated composite skew shells
using higher order shear deformation theory, Lat. Am. J. Solids Struct. 10 (5)
creases as the skew angle, the aspect ratio and the hub radius (2013) 891–919.
ratio of 𝜎 > 1 increases, while it increases as the hub radius ratio [19] A. Kumar, P. Bhargava, A. Chakrabarti, Vibration of laminated composite skew
of 𝜎 < 1 increases. The buckling behaviors will be eliminated hypar shells using higher order theory, Thin-Walled Struct. 63 (63) (2013) 82–90.
with the situation of a non-zero stagger angle. [20] R. Thendiyath, Buckling of laminated composite skew plate using FEM and
machine learning methods, Eng. Comput. 38 (1) (2021) 501–528.
[21] A.A. Shabana, Definition of the slopes and the finite element absolute nodal
CRediT authorship contribution statement coordinate formulation, Multibody Syst. Dyn. 1 (3) (1997) 339–348.
[22] A.M. Mikkola, A.A. Shabana, A non-incremental finite element procedure for
Jing Zhang: Writing – original draft, Software, Data curation, Con- the analysis of large deformation of plates and shells in mechanical system
applications, Multibody Syst. Dyn. 9 (2003) 283–309.
ceptualization. Xiaokang Du: Writing – review & editing, Validation.
[23] O.N. Dmitrochenko, D.Y. Pogorelov, Generalization of plate finite elements for
Yuanzhao Chen: Funding acquisition. Xian Guo: Writing – review absolute nodal coordinate formulation, Multibody Syst. Dyn. 10 (1) (2003)
& editing, Funding acquisition. Liang Li: Writing – review & editing, 17–43.
Funding acquisition. Dingguo Zhang: Writing – review & editing, [24] K. Dufva, A.A. Shabana, Analysis of thin plate structures using the absolute nodal
Validation, Project administration, Funding acquisition. coordinate formulation, Proc. Inst. Mech. Eng. K 219 (4) (2005) 345–355.
[25] C. Liu, Q. Tian, H.Y. Hu, New spatial curved beam and cylindrical shell elements
of gradient-deficient absolute nodal coordinate formulation, Nonlinear Dynam.
Declaration of competing interest 70 (3) (2012) 1903–1918.
[26] D. Yan, C. Liu, Q. Tian, et al., A new curved gradient deficient shell element
of absolute nodal coordinate formulation for modeling thin shell structures,
The authors declare that they have no known competing finan-
Nonlinear Dynam. 74 (1) (2013) 153–164.
cial interests or personal relationships that could have appeared to [27] A.I. Valkeapää, H. Yamashita, P. Jayakumar, et al., On the use of elastic middle
influence the work reported in this paper. surface approach in the large deformation analysis of moderately thick shell
structures using absolute nodal coordinate formulation, Nonlinear Dynam. 80
(3) (2015) 1133–1146.
Data availability
[28] A.A. Shabana, C.J. Desai, E. Grossi, et al., Generalization of the strain-split
method and evaluation of the nonlinear ANCF finite elements, Acta Mech. 231
Data will be made available on request. (4) (2020) 1365–1376.

15
J. Zhang, X. Du, Y. Chen et al. Thin-Walled Structures 188 (2023) 110840

[29] Q.T. Wang, Z. Zhao, Y. Zhang, et al., Nonlinear large deformation of graphene [34] H.H. Yoo, C. Pierre, Modal characteristic of a rotating rectangular cantilever
sheets using a nonlocal-based ANCF method, Int. J. Appl. Mech. 11 (04) (2019) plate, J. Sound Vib. 259 (1) (2003) 81–96.
1950037. [35] K.M. Liew, Y. Xiang, S. Kitipornchai, et al., Vibration of thick skew plates
[30] P. Lan, Y.Q. Cui, Z.Q. Yu, A novel absolute nodal coordinate formulation thin based on mindlin shear deformation plate theory, J. Sound Vib. 168 (1) (1993)
plate tire model with fractional derivative viscosity and surface integral-based 39–69.
contact algorithm, Proc. Inst. Mech. Eng. K 233 (3) (2019) 583–597. [36] W.H. Liu, W.C. Chen, Vibration analysis of skew cantilever plates with stiffeners,
[31] Y. Feng, J.M. Wu, M.Y. Shao, Dynamic characteristics analysis of high speed J. Sound Vib. 159 (1) (1992) 1–11.
rotating thin film antenna, Results Phys. (2022) 105612.
[37] L. Zhou, W.X. Zheng, Vibration of skew plates by the MLS–Ritz method, Int. J.
[32] Y.Z. Chen, D.G. Zhang, L. Li, Dynamics analysis of a rotating plate with a setting
Mech. Sci. 50 (7) (2008) 1133–1141.
angle by using the absolute nodal coordinate formulation, Eur. J. Mech. A Solids
[38] O.G. McGee III, Flexural vibrations of clamped-free rhombic plates with corner
74 (2019) 257–271.
stress singularities. Part II: comparison of solutions, J. Vib. Control 21 (13)
[33] J. Zhao, Q. Tian, H. Hu, Modal analysis of a rotating thin plate via absolute
nodal coordinate formulation, J. Comput. Nonlinear Dyn. 6 (4) (2011) 041013. (2015) 2639–2660.

16