Closed-loop subsonic flutter control of smart variable stiffness composite

ed
plates in tapered-swept configurations

Pritam Mondal, Prashanta K. Mahato*

Department of Mechanical Engineering, IIT (ISM) Dhanbad, Jharkhand 826004, India

iew
Abstract

v
This study investigates the free vibration, transient dynamic response, and flutter
characteristics of cantilevered smart variable-stiffness composite laminated (VSCL) plates,

re
focusing on both forward-swept and backward-swept tapered configurations with identical
aspect ratios. A key objective is implementing closed-loop subsonic flutter control using active
fiber composites (AFCs) as actuators and sensors, incorporating a full-state feedback
er
controller. Additionally, the effect of the plate’s geometry on the flutter characteristics of
VSCL plates is examined. A MATLAB-based finite element (FE) code is employed for normal
pe
mode analysis. Unsteady aerodynamic forces are calculated using the doublet lattice method
(DLM) in MSC Nastran for flutter analysis under subsonic flow conditions. The aerodynamic
model is integrated with structural modal parameters via direct matrix abstraction programming
(DMAP) in Nastran. A simple and effective methodology is used for meshing swept tapered
ot

plates. Aerodynamic loads are interpolated for flutter conditions using MATLAB, and the
controller is designed and simulated in Simulink for closed-loop vibration/flutter control.
Results show that forward-swept plates are prone to divergence, while backward-swept plates
tn

exhibit controlled flutter behavior with the implemented controller.

rin

Keywords: Variable stiffness composite laminate (VSCL); Swept plates; Active fibre
composite (AFC) layer; Pole placement controller; Transient dynamic analysis; Subsonic
ep

flutter
Pr

*Corresponding author.
E-mail addresses: pritamm32@[Link] (P. Mondal), pkmahato@[Link] (P. K. Mahato)

This preprint research paper has not been peer reviewed. Electronic copy available at: [Link]
 1. Introduction

ed
Composite materials are widely used in aerospace, mechanical, and civil engineering
sectors due to their advantages, including low weight, corrosion resistance, and high fatigue
strength [1–3]. Their exceptional strength-to-weight ratio and superior fatigue resistance allow

iew
them to outperform traditional materials in dynamic, weight-sensitive applications. A
significant advancement in this area is the development of Variable Stiffness Composite
Laminates (VSCLs), which employ curvilinear fiber paths to adjust stiffness according to
specific requirements. This innovation has the potential to create structures that are both more

v
efficient and stronger than conventional constant stiffness composite laminates (CSCLs) with
straight fibers [4].

re
Plates are fundamental structural components used across various engineering
applications, such as aircraft wings, fuselages, vehicle bodies, and building slabs, where
er
efficient material usage and structural integrity are critical. Their analysis often serves as an
essential first step before progressing to more complex geometries in structural design.
Research has extensively examined straight plates, both variable and constant stiffness,
pe
particularly regarding their free vibration behavior [5–8]. Understanding these dynamic
characteristics is crucial for analyzing transient responses, where dynamic loads can
significantly impact structural integrity and overall plate behavior. Piezoelectric actuators and
sensors are vital in controlling these transient responses [9–11], providing effective vibration
ot

control solutions. Dynamic instability in structures can also arise from aerodynamic loads, with
flutter being a prominent phenomenon. Flutter occurs when aerodynamic forces interact with
tn

an object's natural structural vibrations, resulting in self-excited oscillations that can amplify
over time. If not addressed, flutter can lead to significant structural damage or failure. Extensive
research has focused on understanding this phenomenon, leading to the development of
rin

effective control strategies and suppression techniques for both subsonic [12–16] and
supersonic [17–19] flow regimes. Collectively, these studies enhance the overall performance
and stability of straight plates under various operating conditions.
ep

While significant advancements have been achieved in understanding the dynamic

behavior and control of straight plates, many engineering applications require the use of non-
rectangular geometries, such as swept, trapezoidal, or skew plates. These designs are
Pr

particularly prevalent in advanced structures like aircraft wings and bridges, where they often
meet specific performance requirements that straight plates may not fulfill. Therefore, studying

This preprint research paper has not been peer reviewed. Electronic copy available at: [Link]
 the behavior of these non-rectangular plates is essential. Early research primarily focused on
isotropic and constant stiffness composite laminates (CSCLs), laying the groundwork for future

ed
advancements in the field, including studies on variable stiffness composite laminates
(VSCLs). Orris and Petyt [20] made notable contributions by implementing high-precision
conforming plate bending elements—both quadrilateral and triangular—within the finite

iew
element method to examine the free vibration characteristics of isotropic triangular and
trapezoidal plates. Subsequently, Srinivasan and Babu [21] conducted a numerical study on the
free vibration of isotropic cantilever quadrilateral plates of general shape. They extended their
work in ref. [22] to analyze flutter using piston theory in isotropic quadrilateral plates under

v
various boundary conditions. Later, in ref. [23], they investigated both free vibration and

re
supersonic flutter of laminated composite quadrilateral plates with clamped edges. In all their
studies, the differential equations of motion were derived in quadrilateral coordinates and
solved using an integral equation technique. The focus on skewed isotropic and laminated
er
structures was further advanced by Yang and Wu [24], who developed a 48-degree-of-freedom
skewed quadrilateral thin shell finite element for static and dynamic analyses of isotropic thin
shell structures, incorporating geometric non-linearity. Utilizing tensorial mathematics, they
pe
modeled a variety of general thin shell structures with arbitrary shapes, including square and
rhombic plates, cylindrical and spherical shells, as well as trapezoidal flat and curved plates.
Their extensive validation through numerous examples demonstrated that skewed meshes
perform comparably to traditional rectangular meshes in both linear and nonlinear static
ot

analyses and in free vibration analyses of plates and shells. Building on this foundation,
Pidaparti and Chang [25] applied this finite element method to account for aerodynamic effects
tn

in supersonic flow, enabling a thorough analysis of flutter in both isotropic and laminated
composite panels. Their study systematically examined the influence of various factors on
flutter boundaries, including boundary conditions, flow angles, fiber orientations, and the
rin

presence of cracks. In a separate exploration, Chowdary et al. [26] used the finite element
method to investigate the dynamic instability of laminated composite skew panels subjected to
supersonic flow. A notable aspect of this study was the modification of the conventional finite
ep

element formulation for rectangular planforms by imposing skew boundary conditions,

allowing for effective modeling of skew panels using isoparametric elements. Barik and
Mukhopadhyay [27] presented a novel four-noded plate bending element developed based on
Pr

Bernoulli's assumptions, addressing shear locking issues and numerical difficulties

encountered by isoparametric elements in analyzing very thin isotropic plates. The proposed
element effectively models arbitrary shapes and simplifies the analysis of free vibration
3

This preprint research paper has not been peer reviewed. Electronic copy available at: [Link]
 behavior without requiring special considerations for curved geometries. The mesh division
becomes straightforward by transforming the arbitrary plate shape into a square domain using

ed
cubic serendipity shape functions, avoiding the complexities associated with irregular
boundaries. In ref. [28], Singha and Ganapathi conducted an analysis of the large amplitude
free flexural vibration of thin laminated composite skew plates using four-noded shear flexible

iew
quadrilateral elements. They extended their research (ref. [29]) to investigate the supersonic
panel flutter behavior of laminated composite skew plates, examining the effects of skew angle,
lay-up configurations, and various boundary conditions on flutter performance. Research on
the vibration characteristics of composite plates has expanded with the integration of smart

v
materials for enhanced control. Kanasogi and Ray [30] investigated Active Constrained Layer

re
Damping (ACLD) treatment to manage nonlinear vibrations in smart skew laminated
composites. ACLD combines a constraining piezoelectric layer with a constrained viscoelastic
layer, significantly improving damping capabilities. A simple velocity feedback control law
er
was employed to implement active damping. They developed a finite element model to analyze
various patch configurations and placements, assessing the effectiveness of different patch
shapes, including rectangular and skew patches. A notable conclusion was that skew patches
pe
outperform rectangular patches when positioned at the non-skewed edges of the substrate,
regardless of skew angle and boundary conditions. This research underscores the significant
role of smart structures in advancing vibration control in composite systems.
ot

While much of the earlier work relied on finite element methods and integral equation
techniques for vibration and flutter analyses, alternative approaches, such as the Rayleigh-Ritz
method, have significantly contributed to the understanding of skewed plate structures. Liew
tn

[31] employed the Rayleigh-Ritz method to analyze the vibrational behavior of symmetrically
laminated cantilever trapezoidal plates. By utilizing two-dimensional orthogonal polynomials
as admissible functions, Liew examined the effects of fiber orientation, layer numbers, and
rin

plate aspect ratios on natural frequencies and mode shapes. Separately, Wang [32] introduced
a B-spline Rayleigh-Ritz method for analyzing the free vibrations of thin skew fiber-reinforced
composite laminates under various boundary conditions. Wang's work addressed the
ep

complexities of arbitrary layups and material anisotropy, demonstrating the method's accuracy
and efficiency through numerical applications.
Pr

Transitioning to variable stiffness laminates, Houmat [33] addressed the gap in

nonlinear free vibration analysis by employing a skew p-element based on thin plate theory
and von Karman strains to study variable stiffness skew laminates. Khalafi and Fazilati [34]
4

This preprint research paper has not been peer reviewed. Electronic copy available at: [Link]
 advanced this field by developing an isogeometric finite element method that utilized cubic
NURBS (Non-Uniform Rational B-Splines) to construct the geometry of variable stiffness

ed
laminated skew plates and derive the governing equations. This research specifically focused
on flutter characteristics of curvilinear skew laminates under yawed supersonic airflow. Farsadi
et al. [35] performed modal analysis on variable stiffness composite skewed plates with various

iew
boundary conditions, employing the generalized differential quadrature (GDQ) method and
optimizing fundamental frequencies through the genetic algorithm. Lastly, Rahamanian et al.
[36] investigated nonlinear flutter and limit cycle oscillations (LCO) of laminated composite
plates featuring curvilinear fiber paths in supersonic flow. Utilizing first-order piston theory

v
and Green strain-displacement relations, this study effectively captured the flutter and post-

re
flutter behaviors of tapered and skew plates. The governing equations were solved using the
GDQ method and Newmark’s method to analyze the time response, emphasizing the effects of
fiber paths, taper ratio, and skewness on the aeroelastic behavior of variable stiffness plates.
er
The conducted literature survey reveals a significant lack of research on closed-loop
vibration control of VSCL tapered swept plates using AFC sensors and actuators. Additionally,
to the best of the authors’ knowledge, no studies have been identified on the flutter control of
pe
such plates under subsonic airflow. To address this gap, the present study develops a
MATLAB-based finite element (FE) code tailored for swept plates, enabling the analysis of
both free vibration and transient dynamic loading. For closed-loop vibration control due to
ot

transient dynamic load, a pole-placement controller is designed and implemented in the

Simulink environment. The modal stiffness, mechanical load vector, and sensor sensitivity
matrices are calculated using the MATLAB FE code and imported into Simulink to model the
tn

controller. The control gain and offset matrices are computed to place the closed-loop poles in
a stable region, ensuring effective vibration control. The uncontrolled response can also be
generated using Newmark's time integration scheme within the developed MATLAB FE code.
rin

In the flutter analysis, unsteady aerodynamic loads are calculated using MSC Nastran with the
DLM, which is applicable for subsonic flow regimes. The modal stiffness matrix, including the
effects of curvilinear fiber, is generated in MATLAB and imported into MSC Nastran using
ep

DMAP. The p-k method is then applied to solve the flutter problem in the modal domain. For
a particular case, the condition defined by the reduced frequency and velocity, at which flutter
already occurred, is identified. Since aerodynamic loads are functions of Mach number and
Pr

reduced frequency, the load obtained from MSC Nastran is interpolated to determine the load
at the reduced frequency using a MATLAB code, capturing the flutter behavior. This

This preprint research paper has not been peer reviewed. Electronic copy available at: [Link]
 aerodynamic load is then used in the Simulink environment, along with other structural modal
parameters. To implement the closed-loop controller, the governing equation of the p-k method

ed
is modified to include actuator force terms, which are expressed in state-space form. With this
modified model, the flutter response is successfully controlled, demonstrating the effectiveness
of the pole-placement controller in stabilizing the system.

iew
In summary, the key contributions of this study are as follows:

 A simple and effective approach for meshing swept plates is employed. Using
this scheme, a case where a cantilever CSCL plate with a cranked planform is

v
modelled, and the natural frequencies are validated with MSC Nastran results.

re
 Closed-loop transient vibration control of tapered VSCL swept plates using a
pole-placement controller is demonstrated.
 Subsonic flutter control of backward-swept tapered VSCL plates is addressed,
utilizing a pole-placement controller.

er
The p-k method is employed for a type of pseudo-transient flutter analysis,
where an actuation effect term is added to the governing equation in state-space
pe
form. Controlled flutter behaviour is demonstrated, with results shown in a
midpoint tip displacement versus time plot, at specific reduced frequency and
velocity.
ot

This research fills a critical gap in the field by providing an effective methodology for
controlling the vibration and flutter of swept VSCL plates with AFC-based actuation.
tn

2. Theoretical formulation

This study investigates smart, tapered swept composite plates with ‘𝑛𝑘’ layers,
incorporating both backward- and forward-swept geometries (see Fig. 1). The top layer
rin

functions as an actuator and the bottom layer as a sensor, both using AFCs, while the
intermediate substrate layers consist of composite materials. The displacement fields are
governed by First-Order Shear Deformation Theory (FSDT) and expressed in matrix form in
ep

(Eq. (1)) as,

Pr

This preprint research paper has not been peer reviewed. Electronic copy available at: [Link]
 U 0 
 
U  1 0 0 Z 0   V0  (1)

ed
U   V   0 1 0 0  
Z  W0    Z  u
W  0 0 1 0 0    x 
    
 y 

iew
with mid-plane displacements U0, V0, and W0 in the X, Y, and Z directions, and rotations of the
normal in the YZ and ZX planes denoted by 𝜃𝑥 and 𝜃𝑦 respectively. The plates are cantilevered
i.e., at Y=0, 𝑈0 = 𝑉0 = 𝑊0 = 𝜃𝑋 = 𝜃𝑌 = 0 (Fig. 1). Finite element discretization of the
domain is achieved using 2D eight-noded quadrilateral isoparametric elements.

v
2.1. Strain-displacement relations

re
The linear strain-displacement relations as shown below are incorporated in the
analysis. er
U V V U V W W U
 xx  ,  yy  ,  xy   ,  yz   ,  xz   (2)
X Y X Y Z Y X Z
pe
ot
tn
rin

Fig. 1. Top view of smart VSCL: (a) backward-swept tapered plate, (b) forward-swept tapered plate

2.2. Coordinate Mapping and Mesh Generation for Swept Tapered Plates
ep

This section describes the mesh generation process for a swept, tapered plate, using the
backward-tapered configuration as an example (Fig. 2). Although both forward- and backward-
tapered configurations are analyzed, only the backward case is explained for clarity. Eight-
Pr

noded isoparametric elements are used to discretize the plate. Unlike meshing a straight

This preprint research paper has not been peer reviewed. Electronic copy available at: [Link]
 rectangular plate, which is straightforward, the described scheme simplifies mesh generation
for a non-rectangular, swept, and tapered plate.

ed
v iew
re
er
pe
ot
tn

Fig. 2. Coordinate mapping and mesh generation (a) Mesh in the computational domain (parametric 𝑠′
― 𝑡′coordinates) (b) Mesh mapped to the physical domain (physical X-Y coordinates) (c) Element in
rin

the physical domain (d) Element mapped to isoparametric domain (local ζ-η coordinates)

First, meshing is conducted within a unit square in the computational 𝑠′ ― 𝑡′ coordinate

system, as shown in Fig. 2(a), which is inherently straightforward. Once meshed, the coordinate
ep

points generated in this domain are mapped onto the physical domain X-Y of the plate (Fig.
2(b)) using the transformation outlined in Eq. (3). In Ref. [37], a similar transformation was
employed to convert the physical domain into the computational domain for solving the
Pr

governing equations and boundary conditions numerically using the differential quadrature
method. However, this study uses the transformation to generate the nodal coordinates of swept
plates by mapping the already meshed unit square plate’s coordinates. This approach allows
8

This preprint research paper has not been peer reviewed. Electronic copy available at: [Link]
 for generating nodal coordinates for arbitrary swept and tapered plates simply by specifying
the geometric parameters (𝑎𝑏, 𝑏ℎ, 𝛼,𝛽), making the method flexible for various designs. Such a

ed
transformation has not been previously applied for nodal coordinate generation in finite
element (FE) analysis in this manner.
Backward swept tapered plate

iew
X  sab  bh (sin  - cos  tan  ) st   t bh cos  tan 
(3)
Y  t bh cos 

In this mapping process, points 𝑃′ in the 𝑠′ ― 𝑡′ system corresponds to points P in the

v
physical X-Y system, representing the same location in different coordinate systems. This
transformation enables the division of the plate into suitable elements and establishes Cartesian

re
coordinates for the nodes on the plate boundary. After mapping and mesh generation in the
physical domain, the (X, Y) coordinates of the eight boundary nodes of each element allow us
to map each element onto a standardized isoparametric domain, defined within the
er  -η plane.
Fig. 2 illustrates one such element's transformation from the physical domain (Fig. 2(c)) to the
isoparametric domain (Fig. 2(d)). Here, point 𝑃 represents the same physical point as P but in
pe
isoparametric coordinates. The shape functions for the isoparametric formulation are defined
in Eq. (4), and the location of any point within each element in the X-Y system is related in
terms of nodal coordinates in Eq. (5). Following this, stiffness and mass matrices are computed
using Gaussian quadrature and then assembled to obtain global data.
ot

i 
1 2
4  
 i  2  i ,  (i  1,2,3,4)
1 2

i   i 1   2 ,   (i  5,7) (4)
tn

4
1
 
i  1  2  2  i ,
4  (i  6,8)
rin

8 8
X i Xi , Y  iYi
i 1 i 1
(5)

For the backward-tapered case (Figure 2(b)), four corner node points are labeled as O(0,0), A
ep

(𝑎𝑏,0), B(𝑎𝑏 + 𝑏ℎ𝑠𝑖𝑛𝛽,𝑏ℎ𝑐𝑜𝑠𝛽), C(𝑏ℎ𝑐𝑜𝑠𝛽𝑡𝑎𝑛𝛼,𝑏ℎ𝑐𝑜𝑠𝛽). When switching to the forward-

tapered configuration, -α, -β are used in Eq. (3) to get Eq. (6), which results in a new set of
points O(0,0), A(𝑎𝑏,0), B(𝑎𝑏 ― 𝑏ℎ𝑠𝑖𝑛𝛽,𝑏ℎ𝑐𝑜𝑠𝛽), C ( ―𝑏ℎ𝑐𝑜𝑠𝛽𝑡𝑎𝑛𝛼,𝑏ℎ𝑐𝑜𝑠𝛽), corresponding to
Pr

𝑜′𝐴′𝐵′𝐶′ in the 𝑠′ ― 𝑡′ domain shown in Fig. 2(a).

This preprint research paper has not been peer reviewed. Electronic copy available at: [Link]
 ed
Forward swept tapered plate
X  sab - bh (sin  -cos  tan )st - tbh cos  tan
(6)
Y  tbh cos 

iew
For both configurations, the base width 𝑎𝑏, the vertical height 𝐿ℎ, and the right slant side length
𝑏ℎ are specified. This transformation approach enables mesh generation tailored to arbitrary
forward/backward-tapered swept plates by generalizing the meshing strategy. The method can

v
be easily adapted to various geometries, such as tapered and curved domains, by merely
modifying the mapping equations while retaining a consistent workflow for FE analysis using

re
conventional eight-noded isoparametric quadrilateral elements. However, care should be taken
to avoid excessive distortion or interior corner angles exceeding 180 degrees, as these can affect
accuracy and solution convergence [27]. er
2.3. Reference fiber path for swept variable stiffness lamina
pe
The reference fiber path for swept laminates with curvilinear fibers is defined by the
following equation:
2(T1  T0 ) b (7)
T (Y )  T0  Y  h cos  for 0  Y  bh cos 
bh cos  2
ot
tn
rin
ep
Pr

Fig. 3. Reference path of the curvilinear fiber

10

This preprint research paper has not been peer reviewed. Electronic copy available at: [Link]
 T0 and T1 represent the fibre angles with respect to the Y-axis at the middle (Y=𝐿ℎ/2) and at
the edges (Y=0 and Y=𝐿ℎ) of the plate respectively. < T0/T1 > is used to represent the fibre path

ed
in a variable stiffness lamina (Fig. 3). To avoid manufacturing issues such as fiber kinking, the
permissible fiber curvature (Ref. [33]) is limited to 3.28 m−1.

iew
2.4. Constitutive relations

The stress-strain relation for the kth layer of the variable stiffness plate in the global (X,
Y, Z) coordinate system is given as:

  Q(Y )  (8)

v
k k k

re
The elements of the transformed stiffness matrix ( Q(Y )  ) for the kth layer are defined in the
k

reference [11].
er
The constitutive relations for the AFC in material coordinates are given by

 p11  Q p11 Q p12 0 0 0   p11   e p11 0 0 

 p   p     p 
 22  Q 12 Q 22 0 0 0   p 22  e p12 0 0
p
 E 1 
pe
 p    
 12   0 0 Q p 66 0 0   p12  -  0 0 0  0 
 p      
 23  
0 0 0 Q p 44 0   p 23   0 e p 42 0  0 
 p13   0 0 0 0 Q p 55   p13   0 0 p 
e 53 
T (9)
 e p11 0 0   p11 
ot

 D p1  e p12 0 0


 p   p
 22   11 0 0   E p1 
 p    p    
D 2   0 0 0   12    0  p 22 0  0 
D p   0 e p 
tn

0  p   0 0  p 33   0 
 3  42   23  
 0 0 p 
e 53   p13 

rin

here, {𝐸𝑝}: electric field vector, [𝑒𝑝]: coefficient matrix of piezoelectric stress, [𝜅𝑝]: dielectric
constants matrix, {𝜀𝑝}: strain vector, [𝑄𝑝]: stiffness matrix, {𝜎𝑝}: stress vector. Superscript p-
denotes piezoelectric-related terms.
ep

The electric field vector ({𝐸𝑝}) of the actuator layer is assumed to be acting along the
X-direction and is expressed as follows
T
 1  
E 
Pr

p
  0 0  V0 (10)
 het 

11

This preprint research paper has not been peer reviewed. Electronic copy available at: [Link]
 
where, ℎ𝑒𝑡 is the spacing between the interdigitated electrodes, and V0 is the potential

ed
difference between two adjacent electrodes. The piezoelectric fibers embedded within the
matrix of an AFC are oriented along the X-axis.

2.5. Finite element modelling

iew
Eight-noded isoparametric quadrilateral elements are employed in the FE formulation
of the plate structure, with each node possessing five independent mechanical degrees of
freedom (𝑈0,𝑉0, 𝑊0 , 𝜃𝑋, 𝜃𝑌). Two additional dependent degrees of freedom, 𝑉𝑎 and 𝑉𝑠, are

v
included in the AFC portion. The elemental mechanical and electrical fields are then
isoparametrically interpolated as follows:

re
8 8 8 8 8
U0  iUi , V0  iVi , W0  iWi , X  iXi , Y  iYi
i 1 i 1 i 1 i 1 i 1
(11)
 8   8 
Va  iVai , Vs  iVsi er
i 1 i 1

The elemental displacements can be written in the matrix form as

pe
U 0  U 
   0 0 0 0   0i 
 V0   
0  0 0 0   V0i 
  8   
W0     0 0  0 0  W0i  (12)
  i 1  0 0 0  0   
x     X i 
   0 0 0 0  
  y     Y 
ot

 i 

By substituting Eq. (11) in Eq. (2), the following relations are obtained
tn


 
8 8
     J   Bu  j u j , E p    J   Bp  j V  , for j thelement
1 1
j
i 1 i 1 (13)
  
  
T
V  Va Vs
rin

where, and [ 𝐽 ], is the Jacobian matrix. [𝐵𝑝] and [𝐵𝑢] are the shape function derivative matrices
for electrical and mechanical fields.
ep

2.6. Governing equations

By using Hamilton’s variational principle, the following (see ref. [12]) governing
equations in terms of global variables can be obtained as follows
Pr

12

This preprint research paper has not been peer reviewed. Electronic copy available at: [Link]
  
 M g  u
     K dd  u   K da  Va   K ds  Vs  Fm
            

   

ed
 K  u   K aa  Va  Fp (14)
 ad   

 K  u   K ss  V    0
 sd    s

iew
The above global coefficient matrices can be formed by assembling the elemental matrices of
each element [38,39]. For 𝑗𝑡ℎ element, these matrices can be defined as
1 1
 M gj        I    u  J d  d ,  I   inertia matrix
T
  u
1 1

v
1 1
 K ddj      B   D  B  J d d ,  D   material constitutive matrix
T
  u c e c

re
1 1
1 1
T T
 K daj    K dsj     B  e p   Be  J d  d   K adj    K sdj 
T
    u    
1 1
1 1 (15)
K     B     Be  J d  d   K 
j T p j
 aa e   ss

1 1
1 1
er
      q  J d d ,
Fmj  u
T
st qst   mechanical surface traction vector
1 1

pe
1 1
 V011p 
F      
p
j
e
T

 het 
 J d  d , piezo - actuator force
1 1

Material constitutive matrix,  Dc  is defined as,

ot

  A Y    B Y   0 
    
 Dc     B Y   D Y  
  0 

(16)
tn

 0 0  A Y   
  

where,
rin

nk Zk

 A , B , D     Q Y  1, Z , Z  dZ ,

ji ji ji
k 1 Z k 1
k
ji
2
( j, i  1, 2,6)

Zk
ep

nk
 Aji    5 / 6     Q ji Y  dZ , ( j, i  4,5)
k
 
k 1 Z k 1

In this analysis, the shear correction factor is considered to be 5/6. The stiffness matrices
Pr

depend on Y, and the calculation of fiber orientation, which is crucial for numerical integration

13

This preprint research paper has not been peer reviewed. Electronic copy available at: [Link]
 during the evaluation of the stiffness matrix, is determined using Eq. (7), specifically at the
Gauss quadrature points.

ed
The Eq. (14) can be rewritten as

M g  u
    Kme  u 
   F   F  (17)

iew
 m a

and,  Fa  and  K me  can be defined as follows

1 1
 Kme   Kdd   Kds  Kss   Ksd    Kda   Kaa  
 Kad 

Fa   Kda  Kaa  Fp 

1

v
Fs    Kds {Vs }   Kds  Kss   Ksd u ,
1
force equivalent of sensor voltage

re
The following transformation reduces the dynamic equation to its modal form. {𝑞ℎ} is
the modal displacement vector and [Λ] is the modal matrix containing eigenvectors of the
system.
u     q h 
er (18)
The modal form of the equation of motion is
pe
qh    K me  qh   Fm  Fa  (19)

where,

K me      Kme   , M g     Mg     I 

-1 -1
ot

 
F    F  ,{F }    F 
m
-1
m a
-1
a

F    F   c q  , c   sensor sensitivity matrix

-1
tn

s s s h s

2.7. Control strategy

rin

In this study, a pole placement full-state feedback controller is implemented to stabilize

the system by shifting unstable poles into the stable half-plane. The control law applied here is
given by:

F    Kctrl i2i  Xˆ  r

ep

a i1   N 
2i1 ii i1

where, (20)

 Kctrl i2i   K1c ii  K2c   ,

ii 
Pr

14

This preprint research paper has not been peer reviewed. Electronic copy available at: [Link]
  
where Fa represents the actuator force, and ‘i’ is the number of modes considered. Kctrl 

ed
is the state feedback gain matrix, and  N  is the offset term that ensures the elimination of

steady-state error, allowing the system to track the desired reference and maintain equilibrium,
even in the presence of disturbances, thus ensuring that the actual output matches the desired

iew

output at the steady state. In this application, the reference r is a zero vector, corresponding

to the cantilever plate's desired equilibrium state. Zero transverse deflection represents the
undisturbed position of the plate, indicating that the system is at rest and no control action is

v
required to maintain this state. While designing the pole-placement full-state feedback
controller, it is ensured that the controllability and observability conditions are satisfied.

re
MATLAB’s ‘place’ command is utilized to assign closed-loop poles and compute

values for Kctrl  based on the desired pole locations. In this setup, the poles are placed in the
er
left half of the complex plane (negative half-plane) to ensure stability, with a typical shift of,
for example, 10 or 20 units. While the closed-loop poles in the pole placement method can be
selected based on specifications such as maximum overshoot, settling time, and other
pe
performance criteria, the main focus here is on demonstrating effective vibration suppression,
rather than meeting these specific performance goals.

2.7.1. Transient dynamic load control

ot

This section presents the state-space model used to address dynamic loading. The

 
mechanical (modal) force vector Fm , modal stiffness  K me  , and sensor sensitivity matrix
tn

cs  are calculated using MATLAB FE code and imported into Simulink for controller design.
The state-space model is expressed for the ‘i’ number of modes in Eq. (21). The block diagram
rin

of the controller is shown in Fig. 4.

ep
Pr

15

This preprint research paper has not been peer reviewed. Electronic copy available at: [Link]
 ed
iew
Fig. 4. Block diagram of a state feedback vibration controller for a smart plate under transient dynamic
loading

v
 
The state vector X̂ , state matrix  Â , input matrix  B̂  , and output matrix Ĉ 

re
are defined in Eq. (21).

 Xˆ  2i1
  Aˆ 
2 i 2 i
 Xˆ  2i1
  Bˆ   
2ii
Fm
i1
  Bˆ 
2ii
F 
a i1
(21)
  
Yˆ = Fs  
i×1
= Cˆ 
i×2i
Xˆ
2i×1
er
where,
pe
   
T
Xˆ  qh qh
2i1 i1 i1

  0ii  I ii  T
 ,  Bˆ    0ii  I   , Cˆ    cs  
  i2i  ii  0ii 
 Aˆ  
  2i2i     K  0ii    2ii  ii 
ot

  me  ii

Uncontrolled responses are generated using the in-house MATLAB code that applies
tn

Newmark's time integration method ([40]) to solve for transient modal responses, considering
the first four modes of vibration. Simulink can also generate an uncontrolled response by
removing the control action, enabling a comparison between vibration suppression with and
rin

without the controller. With structural damping neglected, the poles of the system are
marginally stable, with zero real parts. The controller shifts these poles to the left by a specified
amount to ensure stability.
ep

2.7.2. Subsonic flutter control

For flutter control analysis, the external (mechanical) force term Fm is replaced with  
Pr

the aerodynamic force term Faero , as described in Eq. (22). 

16

This preprint research paper has not been peer reviewed. Electronic copy available at: [Link]
 1  a crVair Qh  ii  ˆ
I

 
Fm    a Vair  QhR  ii
2
   
 X   Daero i2i Xˆ   Faero  (22)

ed
 2 4 kred 


This term is then substituted into the system model (Eq. (21)), resulting in an updated state-
space model (Eq. (23)). The updated state matrix,  Aˆ a  , which accounts for the contribution

iew
 

from the aerodynamic load, is presented in this section. Vair , a , cr and k red represent the flow

velocity, air density, reference chord length, and reduced frequency, respectively. QhR  , QhI 

are the real and imaginary parts of the aerodynamic force matrix (in the modal domain). This

v
modification allows us to analyze the behavior of the system under subsonic airflow and to
apply the control strategy effectively.

re
 Xˆ  2i×1
 Xˆ    Bˆ  F 
  Aˆ a 
2i×2i 2i×1 2ii
a i1

(23)
Yˆ = F  = Cˆ   Xˆ 
s i×1
i×2i 2i×1
er
where,
  0 I  
pe
 
 Aˆ a      a Vair  QhR  
2
 a crVair QhI  
      K    
  me  2  4 kred 
   

The block diagram of the full-state feedback system is depicted in Fig. 5.

ot
tn
rin

Fig. 5. Block diagram of a state feedback controller for flutter control of a smart plate
ep

By considering {𝑟} as a zero vector, the following form for the modal flutter analysis
can be derived from Eq. (23)
Pr

   QhI      a Vair   QhR    

2
 c V
    K 2  
2 c a r ai r
     K1     K me  
c
   q h   0 (24)
   4 k red    2  
      

17

This preprint research paper has not been peer reviewed. Electronic copy available at: [Link]
 As mentioned earlier, structural damping is neglected in the present analysis. Additionally, to
evaluate the flutter parameters associated with uncontrolled behaviour, the contributions from

ed
[𝐾𝑐1] and [𝐾𝑐2] are omitted from Eq. (24). The resulting equation is then solved using the p-k
method and represented as
  a Vair  QhR   

iew
  c V QhI    2

   
2 a r air        K me     qh   0 (24a)
  4 kred   2 
    

 Aeroelastic modeling

v
Aeroelastic modelling consists of three primary stages: developing an FE model of the

re
plate structure, constructing the aerodynamic model, and implementing aero-structure
coupling. The first four modes are considered in the analysis. In MSC Nastran, the standard
tools cannot directly model VSCL plates or capture the actuation effects of an AFC layer. To
er
address this, a "dummy" FE structural model of a CSCL plate with an arbitrary layup, such as
90/0/90/0/90, is created in MSC Patran, ensuring that the geometric properties remain
pe
consistent with those of the VSCL plate to be modelled. The AFC layer is modelled as an
additional mechanical layer within this structure. One key consideration is that as the number
of layers increases, it alters the effective stiffness distribution, leading to localized vibrations,
especially at higher frequencies. This makes higher-order modes difficult to ignore in the
ot

analysis beyond a certain number of layers. For aerodynamic modelling and aero-structure
coupling, MSC FlightLoads is utilized to define the aerodynamic response, which is analyzed
at varying Mach numbers and reduced frequencies. Given that the structural and aerodynamic
tn

grids are independent, a surface spline technique is applied to map the structural degrees of
freedom to the aerodynamic degrees of freedom, enabling an integrated aeroelastic analysis.
Further details on load generation and splining techniques can be found in the reference [41].
rin

To account for curvilinear fibers and actuation effects, modal mass and stiffness matrices are
generated using a MATLAB-based FE code tailored for swept smart VSCL plates. These
matrices are then exported to the MSC Nastran solver via the Direct Matrix Abstraction
ep

Programming (DMAP) module, facilitating the calculation of the generalized aerodynamic

load matrix ([𝑄ℎ ]), which is a complex quantity. The unsteady aerodynamic load is computed
using the Doublet-Lattice Method (DLM), which approximates the lifting surface as a flat plate,
Pr

neglecting thickness effects. These loads are calculated for a range of 𝑘𝑟𝑒𝑑 and Mach numbers
under sea-level atmospheric conditions (𝜌𝑎=1.225 kg/𝑚3) at a low subsonic velocity (M = 0.5).
18

This preprint research paper has not been peer reviewed. Electronic copy available at: [Link]
 The airflow is parallel to the X -direction (Fig. 1). Since the generalized aerodynamic loads are
discrete, a least-squares curve-fitting technique approximates the load coefficients as

ed
continuous functions over the reduced frequency (𝑘𝑟𝑒𝑑) scale for flutter control applications.
Flutter analysis is conducted in incremental velocity steps using the p-k method, which relies
on interpolated and extrapolated unsteady aerodynamic loads. Finally, MSC Nastran is used to

iew
solve the flutter analysis yielding the flutter velocity and frequency for a comprehensive
assessment of the aeroelastic behavior of the structure.

By solving Eq. (24a), the eigenvalue λ is obtained, which is generally a complex

v
quantity (𝜆 = 𝑅𝑒(𝜆) +𝑖.𝐼𝑚(𝜆)). Re(λ) is associated with damping, 𝑔𝑑=2Re(λ)/Im(λ) and
indicates whether deformation grows or decays, while the imaginary part corresponds to the

re
oscillatory nature, representing the frequency of vibration, 𝑓𝑑 = Im(λ)/(2π) in Hz. A positive
Re(λ) implies oscillation growth (instability), whereas a negative Re(λ) signifies decay
(stability). When the Re(λ) becomes zero while the Im(λ) remains non-zero, it indicates the
er
onset of flutter. Before the onset of flutter, the negative Re(λ) denotes stability; after the onset,
the Re(λ) turns positive, indicating instability. The onset of divergence occurs when the Re(λ)
pe
transitions from zero to positive, marking the boundary between stable and unstable
equilibrium. A negative Re(λ) before this transition signifies stability, while a positive Re(λ)
after it denotes instability. Divergence, being a static instability, is characterized by a zero
Im(λ), as oscillatory behavior is absent. For real roots, the damping, 𝑔𝑑 in this case is expressed
ot

as

cr
gd 
tn

(25)
 ln 2Vair

For the analysis at hand, only the first four modes are considered. When using the
rin

DMAP module to extract aerodynamic loads ([𝑄𝑅h ] and [𝑄𝐼h]) from MSC Nastran, the results
are generally provided for all user-specified Mach numbers and reduced frequencies. By
solving Eq. (24a) iteratively in incremental velocity steps, the system's stability can be assessed
ep

by examining the real and imaginary parts of the eigenvalues. If flutter occurs, the aerodynamic
Im(λ)
loads at that 𝑉𝑎𝑖𝑟 and 𝑘𝑟𝑒𝑑( = 𝑐𝑟. 2.𝑉𝑎𝑖𝑟 ) become crucial parameters for defining the flutter
Pr

condition. Consequently, the matrices, [𝑄𝑅h ] and [𝑄𝐼h] are obtained for the specific 𝑉𝑎𝑖𝑟 and
𝑘𝑟𝑒𝑑 values (associated with flutter) by interpolation or extrapolation using a developed

19

This preprint research paper has not been peer reviewed. Electronic copy available at: [Link]
 MATLAB code. Those matrices are used in Eq. (23) to create the state space model in
Simulink.

ed
 Implementation of the controller strategy

With the above in mind, the next steps in the analysis involve designing the control

iew
strategy and assessing the flutter behavior of the system. From the MATLAB FE code, [𝐾𝑚𝑒],
modal structural stiffness matrix (corresponding to the first four modes), and sensor sensitivity
matrix ([𝑐𝑠]) are obtained, and a state-space model is constructed using Eq. (23). The system
and controller are then modeled in the Simulink environment. While designing the pole-

v
placement full-state feedback controller, it is ensured that the controllability and observability

re
conditions are satisfied. The selection of the first four modes is deliberate, as it allows us to
retain the essential system dynamics while ensuring that the reduced-order model remains both
controllable and observable—a critical requirement for effective state feedback control. The
er
state feedback controller with pole placement was selected for its simplicity. Higher modes
were excluded to preserve controllability and observability, as their inclusion would require
more advanced strategies. The focus is on the structural modeling, meshing, and aeroelastic
pe
analysis of VSCL swept tapered plates, with control implemented solely to demonstrate the
controlled response. In the case of flutter, the closed-loop system is designed such that the
open-loop poles are shifted by a certain amount (for example, 20 units) to the left half of the
complex plane to ensure stability. The controller parameters [𝐾𝑐𝑡𝑟𝑙] and [𝑁] are calculated to
ot

ensure that the controller performs as desired, achieving proper closed-loop dynamics and
eliminating steady-state error in reference tracking. The flowchart of the process is shown in
tn

Fig. 6. While poles can be set based on criteria such as settling time and peak overshoot, this
pole-shifting strategy is used for simplicity and to demonstrate the control action. The Simulink
model then generates a plot of transverse deflection versus time, with or without the controller
rin

action, illustrating the system's dynamic response. This is not a full transient analysis, but rather
a kind of pseudo-transient analysis, which serves as a representation to show whether the
response is diverging or not at that condition.
ep
Pr

20

This preprint research paper has not been peer reviewed. Electronic copy available at: [Link]
 ed
v iew
re
er
pe
ot
tn
rin

Fig. 6. Flowchart of the flutter Control Strategy (considering first four modes): MATLAB FE Code,
MSC NASTRAN for aero load computation, and Simulink for control design
ep
Pr

21

This preprint research paper has not been peer reviewed. Electronic copy available at: [Link]
 3. Results and discussions

ed
This section provides an in-depth analysis of free vibration behavior, transient dynamic
vibration control, and subsonic flutter suppression in swept laminates. The accuracy of the
proposed numerical model and approach is first validated through rigorous verification studies.

iew
Following this, the methodology is employed to meet the primary objectives of this study.

3.1. Numerical verifications

A thorough verification of the current numerical model by comparing free vibration

v
and transient response results with literature and results generated using the MSC Nastran

re
solver is presented here. The transient response analysis includes separate cases: one for
transient dynamic loads and another for unsteady aerodynamic loads due to subsonic airflow.

3.1.1. Free vibration analysis er

In this subsection, the free vibration frequencies obtained from the current MATLAB
code are validated through comparisons with established results from multiple sources. The
pe
validation process encompasses three key cases: first, the natural frequencies of skew CSCL
plates are compared with reference values provided by [28], ensuring accuracy in capturing the
effects of plate geometry on vibrational characteristics. Next, the free vibration frequencies of
the CSCL cranked plate are verified against results generated by the MSC Nastran solver,
ot

serving as a benchmark for the present model’s capability to handle complex plate geometries.
Finally, the vibrational frequencies for skew VSCLs are validated against reference [33],
tn

confirming the model's accuracy in simulating composite swept plates with curvilinear fibers.
Each case will be discussed in detail in the subsequent sections.

 CSCL skew plate

rin

A skew-clamped CSCL plate (a=b=1, a/h=1000) of lamination sequence (45°/-

45°/45°/-45°/45°) is considered for the free vibration analysis. The problem descriptions, which
ep

include the configuration and conditions of the skew-clamped plate, along with material
properties, are sourced from the reference [28]. It is observed from Table 1 that the present data
are in good agreement with the reference data. This indicates that the present MATLAB code
Pr

can accurately capture the effects of skew angles.

22

This preprint research paper has not been peer reviewed. Electronic copy available at: [Link]
 Table 1
Validation of non-dimensional natural frequencies of a skew laminated plate

ed
B.C Skew Modes
angle 1 2 3 4 5 6
0° Present 3.9074 7.1586 8.4984 11.2572 13.3676 14.8828
Ref. [28] 3.8996 7.1426 8.4532 11.2008 13.3040 14.7223

iew
CC 30° Present 4.5477 8.3985 9.9074 12.8978 15.7814 17.5593
Ref. [28] 4.5425 8.3787 9.8764 12.8428 15.6673 17.4547
45° Present 6.3163 10.8399 14.5637 15.5287 21.2111 22.2018
Ref. [28] 6.2903 10.7941 14.4487 15.4235 20.9774 21.9865

v
re
 CSCL cranked plate

To validate the proposed coordinate mapping and transformation technique described

in section 2.2 for modeling complex plate geometries, an FE model of a cranked CSCL plate
er
was created in MATLAB, as shown in Fig. 7. The cranked arrow geometry consists of two
backward-tapered, swept quadrilaterals, as defined by Eq. (3). Angles 𝛼1 = 15.243°, 𝛽1= -5°,
pe
𝛼2 = 19.1726°, and 𝛽2 = 10°. The edge OA is fixed. Material properties used are 𝐸1= 128 GPa,
𝐸2= 9.4 GPa, µ12= 0.3, ρ = 1600 kg/𝑚3, 𝐺12 = 𝐺13 = 6.3 GPa, and 𝐺23= 3.15 GPa. Unless
otherwise mentioned, these properties were consistently applied throughout the paper for the
simulations and results to address the current objectives. The plate comprises 8 layers with a
ot

total thickness of 0.006 m. The mapping and transformation scheme was applied in
MATLAB to generate the mesh. Each quadrilateral was mapped as a unit rectangle and then
tn

meshed with 8-noded isoparametric elements. The nodal coordinates were then transformed
from the computational domain back to the physical domain using Eq. (3).
rin
ep
Pr

23

This preprint research paper has not been peer reviewed. Electronic copy available at: [Link]
 ed
v iew
re
er
Fig. 7. Geometry of cranked CSCL plate
pe
Then isoparametric formulation was incorporated to create the FE model of the geometry.
Natural frequencies generated in MATLAB were verified against results from MSC Nastran.
MSC Patran served as the pre-and post-processing platform, while MSC Nastran was used to
ot

solve the problem. To construct the geometry, two quadrilaterals were created using the six
key points, and two surfaces were modeled. Nodes along the common edge (connecting Points
tn

E and B) were merged to ensure continuity, treating them as a single node. The first four mode
frequencies, shown in Table 2, demonstrate that the MATLAB-based results closely align with
those from MSC Nastran.
rin

Table 2
Validation of natural frequencies (Hz) of CSCL cranked plate
Mode 1 Mode 2 Mode 3 Mode 4
(0°/90°)4
ep

MATLAB 19.67 75.83 121.66 210.28

MSC Nastran 19.29 75.73 121.26 207.79
[(0°/90°)2]S
MATLAB 16.77 70.66 104.54 194.01
Pr

MSC Nastran 16.53 70.44 104.36 191.91

24

This preprint research paper has not been peer reviewed. Electronic copy available at: [Link]
  VSCL skew plate

ed
The MATLAB code developed in this study was utilized to compute the natural
frequencies of skew VSCL plates, with the results compared and validated against data
available in the literature. For the free vibration analysis, a clamped square plate (0.5 m × 0.5

iew
m × 0.005 m) featuring a lamination sequence of (<45°/T1>, −<45°/T1>, <45°/T1>) was
analyzed. The material properties were adopted as specified in ref. [33]. In the current
implementation, the domain was discretized into 100 elements. As presented in Table 3, the
calculated frequencies show strong agreement with the reference data, verifying the model’s

v
capability to accurately simulate composite plates with curvilinear fibers.

re
Table 3
Validation of frequency parameter of clamped VSCL skew plate
𝛼(°)
𝑇1(°) 15 30 45

45
Present
.3297
erRef.[33]
.3332
Present
.4421
Ref.[33]
.4496
Present
.6697
Ref.[33]
.6913
65 .3218 .3252 .4269 .4337 .6190 .6344
pe
90 .3198 .3230 .4072 .4128 .5566 .5677

3.1.2. Transient response analysis

ot

To validate the accuracy of our model, two types of cantilever plates were analyzed:
one with a forward sweep (denoted as "frwd") and another with a backward sweep (denoted
tn

as "bkwd"), as illustrated in Fig. 8. The edge OA is fixed. Each CSCL consists of five layers
with a (90°/0°/90°/0°/90°) stacking sequence, resulting in a total thickness of 0.006 m. The
material properties are consistent with those used for the cranked arrow plate case. A downward
rin

sinusoidal force of 100 N was applied at the mid-tip of the cantilever plate: at Mf for the frwd
and at Mb for bkwd (Fig. 8). The force was applied for 0.1 seconds at a frequency of 50 Hz and
then removed, with the total simulation running for 0.5 seconds. Structural damping was
ep

neglected.
Pr

25

This preprint research paper has not been peer reviewed. Electronic copy available at: [Link]
 ed
v iew
re
Fig. 8. Geometry of frwd (O-A-D-E) and bkwd (O-A-B-C) swept plates

er
For the modal transient analysis, the first four modes were considered. The transverse
displacement at the midpoint tip over time, generated using both MATLAB code and MSC
Nastran, is presented in Fig. 9, with the first four natural frequencies summarized in Table 4.
pe
Both the table and the plot include data for bkwd and frwd, which exhibited the same natural
frequencies and transient responses. This similarity can be attributed to the material symmetry
and geometric configuration of the plates. Notably, the choice of 0° and 90° fiber orientations
in the stacking sequence played a crucial role in ensuring dynamic equivalence between bkwd
ot

and frwd. The MATLAB-generated results and the MSC Nastran results showed close
agreement, confirming the accuracy of the modeled plates in capturing the expected dynamic
tn

response.
Table 4
Validation of natural frequencies of frwd/bkwd swept (α = β = 30°) plates
Mode 1 Mode 2 Mode 3 Mode 4
rin

MATLAB 15.84 64.77 124.89 189.39

MSC Nastran 15.82 64.41 123.78 188.08
ep
Pr

26

This preprint research paper has not been peer reviewed. Electronic copy available at: [Link]
 ed
v iew
re
Fig. 9. Transient dynamic load response for frwd/bkwd swept plate
er
Now, the dynamic behavior of frwd and bkwd under subsonic airflow is analyzed, with
pe
flutter characteristics for the first two modes summarized in Table 5 and Table 6. By solving
the equations of motion iteratively with increasing velocity steps, according to Eq. (24a), the
data in these tables were obtained.
ot

For bkwd (Table 5), at Mode 2, there is a transition in damping between 175 m/s and
200 m/s, where the damping crosses zero while the frequency remains non-zero. To determine
the exact flutter velocity, damping, and frequency are plotted against velocity; the flutter point
tn

is identified where damping reaches zero at a non-zero frequency. While such data
visualization could provide additional insights, this study focuses on the time response of
transverse displacement at the midpoint tip of the plate. System dynamics were modeled in
rin

Simulink without the actuator sensor effect using the state-space form in Eq. (23), considering
the first four modes. Displacement data obtained from Simulink is converted to physical
coordinates from modal coordinates to represent a time response plot.
ep
Pr

27

This preprint research paper has not been peer reviewed. Electronic copy available at: [Link]
 Table 5
Flutter summary of the bkwd plate

ed
Reduced Velocity, Damping, Frequency, Eigenvalue (λ)
Mode frequency 𝑉𝑎𝑖𝑟 𝑔𝑑 𝑓𝑑 Real part, Imaginary
(𝑘𝑟𝑒𝑑) (m/s) (Hz) Re(λ) part, Im(λ)

iew
1 0.3056 50 -0.0834 16.2136 -4.2487 101.8731
0.1689 100 -0.1628 17.9162 -9.1635 112.5710
0.1332 150 -0.2447 21.1996 -16.2949 133.2013
0.1281 175 -0.2947 23.7905 -22.0283 149.4801
0.1298 200 -0.3749 27.5446 -32.4456 173.0680

v
0.1261 250 -0.9972 33.4456 -104.7732 210.1446

re
0.1007 300 -1.6768 32.0539 -168.8501 201.4003
2 1.1978 50 -0.0201 63.5428 -4.0146 399.2509
0.5775 100 -0.0332 61.2711 -6.3882 384.9776
0.3585 150 -0.0390 57.0601 -6.9911 358.5193
0.29 175
er
-0.0310 53.8424 -5.2484 338.3018
0.2325 200 0.0077 49.3445 1.1893 310.0405
pe
0.1577 250 0.4947 41.8437 65.0306 262.9115
0.1286 300 0.9333 40.9352 120.0221 257.2031

When 𝑘𝑟𝑒𝑑=0.29 and 𝑉𝑎𝑖𝑟=175 m/s (Table 5), the Re(λ) is negative, causing any initial
displacement to decay as the midpoint tip displacement stabilizes over time (Fig. 10). However,
ot

at 𝑘𝑟𝑒𝑑=0.2325 and 𝑉𝑎𝑖𝑟=200 m/s, the Re(λ) becomes positive, indicating the diverging
oscillatory behavior of flutter, as shown in Fig. 11. These figures also demonstrate that the time
tn

response generated in Simulink effectively captures the physics of the problem, aligning with
the data in Table 5, thereby validating the methodology. For the frwd case (Table 6),
divergence (static instability) occurs between 125 m/s and 150 m/s.
rin
ep
Pr

28

This preprint research paper has not been peer reviewed. Electronic copy available at: [Link]
 ed
v iew
re
Fig. 10. Converging mid-point tip displacement (transverse) vs. time for bkwd plate under aero load at
𝑘𝑟𝑒𝑑=0.29 and 𝑣𝑎𝑖𝑟=175 m/s er
While the natural frequencies and transient dynamic responses under point load are
similar for both geometries, the two configurations exhibit distinct behavior under aerodynamic
pe
load, as expected: the bkwd configuration experiences flutter, whereas the frwd configuration,
which is known to be more prone to divergence, shows divergence at a considerably lower
speed. Due to the non-oscillatory behavior associated with the divergence of the frwd
configuration, only the bkwd’s time response under the flutter condition is presented. This
ot

concludes our validation study, confirming the accuracy of our modeling approach.
tn
rin
ep
Pr

29

This preprint research paper has not been peer reviewed. Electronic copy available at: [Link]
 ed
v iew
re
Fig. 11. Diverging mid-point tip displacement (transverse) vs. time for bkwd plate under aero load at
𝑘𝑟𝑒𝑑=.2325 and 𝑣𝑎𝑖𝑟=200 m/s, indicating flutter
er
Table 6
Flutter summary of the frwd plate
pe
Reduced Velocity, Damping, Frequency, Eigenvalue (λ)
Mode frequency 𝑉𝑎𝑖𝑟 𝑔𝑑 𝑓𝑑 Real part, Imaginary
(𝑘𝑟𝑒𝑑) (m/s) (Hz) Re(λ) part, Im(λ)

1 0.2753 50 -0.0408 14.6047 -1.8717 91.7641

ot

0.1010 100 -0.1192 10.7213 -4.0139 67.3641

0.0460 125 -0.2789 6.1036 -5.3481 38.3501
0.0000 150 0.1251 0.0000 43.3499 0.0000
tn

0.0000 200 0.2302 0.0000 106.3673 0.0000

0.0000 250 0.2730 0.0000 157.6696 0.0000
0.0000 300 0.3004 0.0000 208.2121 0.0000
rin

2 1.2065 50 -0.0279 64.0048 -5.6090 402.1541

0.5914 100 -0.0491 62.7512 -9.6697 394.2773
0.4662 125 -0.0601 61.8283 -11.6731 388.4784
0.3814 150 -0.0711 60.7070 -13.5622 381.4335
ep

0.2728 200 -0.0918 57.8951 -16.7052 363.7656

0.2052 250 -0.1080 54.4243 -18.4644 341.9582
0.1587 300 -0.1157 50.5011 -18.3620 317.3079
Pr

30

This preprint research paper has not been peer reviewed. Electronic copy available at: [Link]
 3.2. Present results

ed
The numerical results related to the study objectives are presented here, organized as
follows: Section 3.2.1 addresses vibration control under transient dynamic loading, and Section
3.2.2 examines the suppression of flutter under subsonic airflow conditions for two forward-

iew
swept (frwd1, frwd2) and two backward-swept (bkwd1, bkwd2) tapered plates, illustrated in
Fig. 12. These plates, with an aspect ratio (A.R) of 2.4, have identical tip and root lengths of
0.2 m and 0.3 m, respectively, and a vertical height of 0.6 m between parallel edges to ensure
a proper comparison of geometric effects on dynamic/aeroelastic behaviour. The

v
corresponding sweep angles, denoted α and β, are provided in Table 7, with plate coordinates
presented in Table 8.

re
er
pe
ot

Fig. 12. Geometries of frwd1 (1-2-8-7), frwd2 (1-2-10-9), bkwd1 (1-2-4-3), and bkwd2 (1-2-6-5)
tn

plates

All four plates share a symmetric laminate stacking sequence of (<90°/60°>, <30°/45°>,
rin

<45°/60°>)s, yielding a total thickness of, ht, 0.0055 m, and two AFC patches, each with a
thickness, hp=0.00025 m, acting as actuator and sensor layers on the laminate's top and bottom
surfaces. Fig. 13 illustrates the computational domain showing the actuator/sensor placement
ep

and corresponding element numbers, using an example where the plate is meshed with 80
elements. The AFC patch, embedded on the surface of the plate, spans 0.5𝑎′𝑏, with a vertical
separation of 0.4𝐿′ℎ between parallel edges.
Pr

31

This preprint research paper has not been peer reviewed. Electronic copy available at: [Link]
 ed
v iew
re
Fig. 13. Computational domain of the plate, highlighting the element numbers comprising the AFC
patch

The electromechanical properties of the AFC material are as follows: the elastic moduli are 𝐸11
er
= 119.7 GPa and 𝐸22 = 129.1 GPa; the shear moduli are 𝐺12 = 39.1 GPa and 𝐺13 = 𝐺23 = 32.4
GPa; Poisson's ratios are µ12 = 0.35 and µ13 = µ23= 0.38; and density is ρ = 6700 Kg/𝑚3. The
pe
dielectric constants are 𝜅11 =8.6 × 10―9 F/m and 𝜅33 = 6.5 × 10―9 F/m, with piezoelectric stress
coefficients 𝑒11 = 14.1 C/𝑚2, 𝑒21 = -3.3 C/𝑚2, and 𝑒24 = 10.8 C/𝑚2. The material properties of
the VSCL are consistent with those used for the cranked CSCL plate. The natural frequencies
ot

for the first four modes are provided in Table 9. The airflow over the plates is subsonic (M=
0.5) and occurs under standard sea-level atmospheric conditions. A detailed analysis of
vibration control due to transient point loading and suppression of flutter under subsonic
tn

airflow is presented in the following subsections.

Table 7
Geometric parameters of frwd1, frwd2, bkwd1, and bkwd2 plates
Configuration ab (m) bh (m) α° β°
rin

bkwd1 .3 .6212 𝛼𝑏1 = 23.49 𝛽𝑏1 = 15

bkwd2 .3 .6620 𝛼𝑏2 = 32.33 𝛽𝑏2 = 25
frwd1 .3 .6542 𝛼𝑓1 = 15 𝛽𝑓1 = 23.49
ep

frwd2 .3 .7101 𝛼𝑓1 = 25 𝛽𝑓1 = 32.33

Pr

32

This preprint research paper has not been peer reviewed. Electronic copy available at: [Link]
 Table 8
Coordinates of frwd1, frwd2, bkwd1, and bkwd2 plates (in the context of Fig. 12)

ed
bkwd1 bkwd2 frwd1 frwd2
Point Coordinates Point Coordinates Point Coordinates Point Coordinates
ID (X, Y) ID (X, Y) ID (X, Y) ID (X, Y)
1 (0,0) 1 (0,0) 1 (0,0) 1 (0,0)

iew
2 (.3,0) 2 (.3,0) 2 (.3,0) 2 (.3,0)
4 (.4608,.6) 6 (.5798,.6) 8 (.0392,.6) 10 (-.0789,.6)
3 (.2608,.6) 5 (.3798,.6) 7 (-.1608,.6) 9 (-.2798,.6)

v
Table 9
Natural frequencies of frwd1, frwd2, bkwd1, and bkwd2 plates

re
Configuration Mode 1 Mode 2 Mode 3 Mode 4
(Hz) (Hz) (Hz) (Hz)
bkwd1 20.99 76.75 119.69 216.19
bkwd2 18.49 75.41
er 114.32 195.89
frwd1 14.10 58.41 104.36 153.19
frwd2 10.98 47.97 100.74 126.11
pe

3.2.1. Control of vibration under transient dynamic load

A downward mechanical force of 100sin(2π × 50t) N is applied at the midpoint of the

ot

tip of four laminated plates. This force is applied for 0.1 seconds and then removed, with the
simulation running for 0.5 seconds. Although the actuators and sensors are located at the plates’
tn

root (Fig. 1)—where the highest bending moments occur—all displacement plots represent the
midpoint tip displacement over time. This method effectively visualizes the control
effectiveness of the overall plate response. Transverse displacement data extracted from the
rin

midpoint tip node is plotted against time using MATLAB finite element code, and the
responses of all four plates are shown in Fig. 14, taking into account four modes for the
transient modal analysis. The displacement values in the physical coordinate system are
ep

converted from the modal displacement data.

Pr

33

This preprint research paper has not been peer reviewed. Electronic copy available at: [Link]
 ed
v iew
re
Fig. 14. Mid-point tip displacement (Transverse) vs. time of frwd1, frwd2, bkwd1, and bkwd2 plates
under transient dynamic load: Uncontrolled response
er
The frwd2 configuration exhibited a second mode frequency of 47.97 Hz (Table 9), which is
close to the forcing frequency of 50 Hz. This proximity resulted in higher amplitude oscillations
and pronounced spikes in the displacement versus time plot. Since structural damping is
pe
neglected, the vibrations do not decay, indicating that all configurations are marginally stable.
To suppress these vibrations, a closed-loop pole placement state feedback controller is
designed. The modal stiffness matrix and modal mechanical force vector generated by the
ot

MATLAB code are incorporated into the state-space model (defined by Eq. (21)) within the
Simulink environment, where the control algorithm is applied (Fig. 4). During the controller
design, observability and controllability conditions are ensured. The reference signal is set to
tn

zero as it represents its equilibrium mean position. The control signal, comprising the controller
gain matrix [𝐾𝑐𝑡𝑟𝑙] and the offset matrix [𝑁], aims to bring the plate back to its mean
equilibrium position. The MATLAB ‘place’ command is employed to position all the closed-
rin

loop poles at desired locations, specifically shifting the marginally stable poles by 10 units into
the left half-plane. This shift ensures that any oscillations will eventually decay over time,
allowing the system to return to its equilibrium position. The controlled response is presented
ep

in Fig. 15.
Pr

34

This preprint research paper has not been peer reviewed. Electronic copy available at: [Link]
 ed
v iew
re
Fig. 15. Mid-point tip displacement (Transverse) vs. time of frwd1, frwd2, bkwd1, and bkwd2 plates
under transient dynamic load: Controlled response
er
3.2.2. Control of flutter at subsonic flow
pe
The section presents the flutter analysis of configurations frwd1, frwd2, bkwd1, and
bkwd2 in subsonic airflow, with the flow direction illustrated in Fig. 1. The methodology
employed aligns with the approach outlined in section 2.7.2. The modal stiffness matrix and
sensor matrix obtained from the MATLAB finite element code, along with the unsteady
ot

aerodynamic load from MSC Nastran, from which [𝑄𝐼ℎℎ]/𝑘𝑟𝑒𝑑 and [𝑄𝑅ℎℎ] at a given reduced
frequency, 𝑘𝑟𝑒𝑑, is obtained by interpolation or extrapolation. These matrices are then fed into
tn

Simulink to create the state-space model (as described by Eq. 23) and to design a pole
placement controller (Fig. 5). The control law (Eq. 20) is applied to return the plate to its set
position, defined as the mean equilibrium position. The analysis considers the first four modes.
rin

Although design variables such as peak overshoot and settling time could be adjusted, this
analysis primarily focuses on demonstrating the controller's effectiveness by shifting all open-
loop poles leftward in the complex plane. This shift ensures that even marginally stable or
ep

unstable poles move into the left-hand plane, leading to decay in any unwanted vibrations. A
flutter summary (of the first two modes) for bkwd1 and bkwd2 is presented in Table 10 and
Table 11, respectively.
Pr

35

This preprint research paper has not been peer reviewed. Electronic copy available at: [Link]
 Table 10
Flutter summary of the bkwd1 plate

ed
Reduced Velocity, Damping, Frequency, Eigenvalue (λ)
Mode frequency 𝑉𝑎𝑖𝑟 𝑔𝑑 𝑓𝑑 Real part, Imaginary
(𝑘𝑟𝑒𝑑) (m/s) (Hz) Re(λ) part, Im(λ)

1 0.3995 50 -0.0729 21.1960 -4.8574 133.1786

iew
0.2127 100 -0.1474 22.5629 -10.4508 141.7669
0.1598 150 -0.2339 25.4329 -18.6920 159.7999
0.1510 200 -0.3901 32.0410 -39.2692 201.3197
0.1543 210 -0.4812 34.3874 -51.9821 216.0626

v
0.1323 250 -1.2478 35.1001 -137.5898 220.5404

re
0.0880 300 -2.4261 28.0041 -213.4392 175.9550
2 1.4188 50 -0.0255 75.2684 -6.0300 472.9251
0.6790 100 -0.0501 72.0463 -11.3346 452.6802
0.4143 150 -0.0711
er 65.9384 -14.7235 414.3034
0.2560 200 -0.0373 54.3347 -6.3719 341.3950
0.2278 210 0.0249 50.7667 3.9696 318.9769
pe
0.1670 250 0.5691 44.3000 79.2084 278.3449
0.1287 300 1.1046 40.9582 142.1290 257.3478

Table 11
ot

Flutter summary of the bkwd2 plate

Reduced Velocity, Damping, Frequency, Eigenvalue (λ)
Mode frequency 𝑉𝑎𝑖𝑟 𝑔𝑑 𝑓𝑑 Real part, Imaginary
tn

(𝑘𝑟𝑒𝑑) (m/s) (Hz) Re(λ) part, Im(λ)

1 0.3546 50 -0.0780 18.8112 -4.6122 118.1944

0.1921 100 -0.1521 20.3847 -9.7378 128.0806
rin

0.1468 150 -0.2244 23.3571 -16.4672 146.7567

0.1356 200 -0.3081 28.7691 -27.8500 180.7617
0.1428 230 -0.4140 34.8572 -45.3383 219.0140
0.1479 240 -0.5169 37.6505 -61.1380 236.5650
ep

0.1485 250 -0.6919 39.3863 -85.6182 247.4711

0.1188 300 -1.5304 37.8017 -181.7485 237.5148
2 1.4010 50 -0.0223 74.3241 -5.2023 466.9922
Pr

0.6844 100 -0.0427 72.6194 -9.7304 456.2814

0.4351 150 -0.0637 69.2552 -13.8667 435.1430

36

This preprint research paper has not been peer reviewed. Electronic copy available at: [Link]
 0.2955 200 -0.0750 62.7175 -14.7754 394.0657
0.2268 230 -0.0300 55.3580 -5.2107 347.8247

ed
0.2043 240 0.0486 52.0352 7.9485 326.9467
0.1875 250 0.1897 49.7305 29.6412 312.4661
0.1479 300 0.7551 47.0729 111.6682 295.7679

iew
For bkwd1, flutter occurs between airspeeds of 200 and 210 m/s, while for bkwd2, it occurs
between 230 and 240 m/s. Notably, both configurations experience flutter in Mode 2. Within
these velocity ranges, the damping coefficient changes sign while maintaining a non-zero

v
frequency, indicating flutter. At 200 m/s for bkwd1 and 230 m/s for bkwd2, the Re(λ) is

re
negative, signifying that initial displacement disturbances will eventually decay. These
responses are not plotted, as the primary goal is flutter control. The similar behavior exhibited
by the bkwd plate is plotted in Fig. 10 of section 3.1.2, showing the time response just before
er
flutter onset. However, at a velocity of 210 m/s with 𝑘𝑟𝑒𝑑 = 0.2278 for bkwd1, and 240 m/s
with 𝑘𝑟𝑒𝑑= 0.2043 for bkwd2, the Re(λ) becomes positive, indicating a diverging amplitude
over time, as illustrated in Fig.16. This indicates flutter, as seen in the transverse midpoint tip
pe
displacement vs. time plot. Fig. 17 presents the controlled response to this flutter behavior.
Table 9 indicates that bkwd1 has a stiffer configuration compared to bkwd2. However, it is
observed that bkwd2, which has larger sweep angles (α, β), exhibits a stabilizing effect with a
ot

higher flutter velocity.

tn
rin
ep
Pr

Fig. 16. Diverging mid-point tip displacement (Transverse) vs. time for bkwd1 and bkwd2 plates under
aero load: Uncontrolled response
37

This preprint research paper has not been peer reviewed. Electronic copy available at: [Link]
 Table 12 and Table 13 present the flutter summaries for frwd1 and frwd2. Both configurations
experience divergence before flutter, as forward-swept plates are more susceptible to this

ed
phenomenon. For frwd1, divergence occurs at Mode 1 between airspeeds of 130 and 140 m/s.

v iew
re
er
Fig. 17. Converging mid-point tip displacement (Transverse) vs. time for bkwd1 and bkwd2 plates
pe
under aero load: Controlled response

At 140 m/s, the reduced frequency 𝑘𝑟𝑒𝑑 is zero, with a positive Re(λ) and zero imaginary part,
indicating non-oscillatory behavior. Under this condition, Eq. (25) is used to compute 𝑔𝑑 in the
ot

flutter summary tables. However, flutter is observed between 250 and 300 m/s at Mode 2. For
frwd2, divergence occurs at a lower speed, between 100 and 110 m/s at Mode 1. The reduced
divergence speed for frwd2, which has larger sweep angles (α, β) than frwd1, indicates
tn

reduced stability against divergence. These results demonstrate the controller's capability to
effectively suppress vibrations across various configurations.
rin

Table 12
Flutter summary of the frwd1 plate
Reduced Velocity, Damping, Frequency, Eigenvalue (λ)
Mode frequency 𝑉𝑎𝑖𝑟 𝑔𝑑 𝑓𝑑 Real part, Imaginary
ep

(𝑘𝑟𝑒𝑑) (m/s) (Hz) Re(λ) part, Im(λ)

1 0.2464 50 -0.1004 13.0697 -4.1239 82.1195

0.0899 100 -0.3125 9.5342 -9.3587 59.9049
Pr

0.0128 130 -2.4894 1.7597 -13.7618 11.0564

0.0000 140 0.0690 0.0000 22.3299 0.0000

38

This preprint research paper has not been peer reviewed. Electronic copy available at: [Link]
 0.0000 150 0.1134 0.0000 39.3002 0.0000
0.0000 200 0.2162 0.0000 99.8899 0.0000

ed
0.0000 250 0.2755 0.0000 159.1176 0.0000
0.0000 300 0.3277 0.0000 227.1360 0.0000
2 1.0841 50 -0.0361 57.5109 -6.5180 361.3513

iew
0.5231 100 -0.0675 55.5035 -11.7681 348.7386
0.3886 130 -0.0848 53.6053 -14.2797 336.8123
0.3558 140 -0.0898 52.8525 -14.9025 332.0822
0.3270 150 -0.0941 52.0409 -15.3768 326.9824

v
0.2226 200 -0.0981 47.2475 -14.5642 296.8645
0.1587 250 -0.0594 42.0976 -7.8610 264.5067

re
0.1188 300 0.0001 37.8124 0.0148 237.5826

Table 13 er
Flutter summary of the frwd2 plate
Reduced Velocity, Damping, Frequency, Eigenvalue (λ)
Mode frequency 𝑉𝑎𝑖𝑟 𝑔𝑑 𝑓𝑑 Real part, Imaginary
pe
(𝑘𝑟𝑒𝑑) (m/s) (Hz) Re(λ) part, Im(λ)

1 0.1815 50 -0.1295 9.6275 -3.9164 60.4913

0.0320 100 -0.8118 3.3957 -8.6598 21.3359
0.0000 110 0.0504 0.0000 12.7999 0.0000
ot

0.0000 150 0.1747 0.0000 60.5343 0.0000

0.0000 200 0.2219 0.0000 102.5240 0.0000
tn

0.0000 250 0.2488 0.0000 143.6953 0.0000

0.0000 300 0.2715 0.0000 188.2151 0.0000
2 0.8927 50 -0.0291 47.3578 -4.3288 297.5577
0.4385 100 -0.0531 46.5249 -7.7650 292.3245
rin

0.3967 110 -0.0579 46.2950 -8.4280 290.8802

0.2837 150 -0.0764 45.1594 -10.8322 283.7447
0.2039 200 -0.0944 43.2666 -12.8306 271.8522
ep

0.1542 250 -0.1004 40.9137 -12.9050 257.0684

0.1200 300 -0.0869 38.1814 -10.4206 239.9011
Pr

39

This preprint research paper has not been peer reviewed. Electronic copy available at: [Link]
 4. Conclusion

ed
This study presents a comprehensive approach for analyzing and implementing closed-
loop control of the tapered swept cantilevered smart VSCL plates with an emphasis on AFC
patches used as both actuators and sensors. The effect of the plate’s geometry on the flutter

iew
characteristics of VSCL plates is also studied. The MATLAB-based finite element (FE)
framework was developed to model free vibration and transient dynamic responses, while MSC
Nastran was used for calculating unsteady aerodynamic loads in subsonic flow conditions.
Additionally, the p-k method was applied to flutter analysis, incorporating an actuation effect

v
term into the governing equation in state-space form. Simulink was employed to design and
implement a closed-loop pole-placement controller for effective vibration/flutter suppression,

re
leveraging the AFC patch for dynamic control. An important contribution of this study is the
generalization of the existing meshing strategy, enabling the modeling of arbitrary swept and
tapered plates. Once the meshing was done, the isoparametric finite element formulation was
er
utilized to calculate the global stiffness matrix, enabling precise structural analysis. This
methodology was used to model a cantilever CSCL plate having cranked planform, with the
pe
results validating against MSC Nastran simulations.

Key findings from this study are:

 Generalized Meshing Approach – Efficiently creates FE models of arbitrary swept and

ot

tapered plates, ensuring geometric flexibility while maintaining accuracy and

computational efficiency. The coordinates of the node points for any arbitrary swept
tn

and tapered plate are generated using a pre-meshed unit square, with the mapping
equation defining the process.

 Effective Transient Vibration Control – The Active Fiber Composite (AFC) patch was
rin

successfully used as both an actuator and a sensor in the closed-loop control system.
The pole-placement controller, coupled with the AFC patch, effectively suppressed
transient vibrations in the VSCL plates, validating its effectiveness in controlling
ep

dynamic responses to transient loads, a crucial capability for practical aerospace

applications.
Pr

 Flutter Behavior and AFC Control – In the forward-swept tapered case, plates with
higher sweep angles (frwd2) showed reduced divergence speeds compared to their

40

This preprint research paper has not been peer reviewed. Electronic copy available at: [Link]
 lower-sweep counterparts (frwd1). The analysis also showed that backward-swept
plates with higher sweep angles (bkwd2) exhibited higher flutter velocities compared

ed
to those with lower sweep angles (bkwd1), indicating the stabilizing effect of increased
sweep. To ensure a proper comparison of geometric effects on aeroelastic behavior, all
plates—both forward- and backward-swept tapered—were designed with identical

iew
aspect ratios, maintaining the same root and tip chord lengths. The controller
demonstrated reliable flutter suppression in backward-swept configurations, indicating
that AFC-based active control offers a promising solution for subsonic flutter control.

v
In summary, this work demonstrated the combined use of MATLAB, MSC Nastran,
and Simulink in modelling, analysing, and controlling the vibration and flutter behaviour of

re
tapered swept composite plates, providing a comprehensive solution for dynamic control in
aerospace applications.
er
CRediT authorship contribution statement
pe
Pritam Mondal: Conceptualization, Writing – review & original draft, Prashanta K.
Mahato: Supervision, Writing – Review & Editing.

References
ot

[1] B. Shue, A. Moreira, G. Flowers, Review of recent developments in composite material for
aerospace applications, In International Design Engineering Technical Conferences and
tn

Computers and Information in Engineering Conference (2009) 811-819,

[Link]
[2] M. Bhong, T.K.H. Khan, K. Devade, B. Vijay Krishna, S. Sura, H.K. Eftikhaar, H. Pal Thethi,
N. Gupta, Review of composite materials and applications, Mater. Today Proc. (2023),
rin

[Link]
[3] V.M. Karbhari, 1-Introduction: the use of composites in civil structural applications, In
Durability of Composites for Civil Structural Applications, Woodhead Publishing Series in
Civil and Structural Engineering (2007) 1–10, [Link]
ep

[4] P. Ribeiro, H. Akhavan, A. Teter, J. Warmiński, A review on the mechanical behaviour of

curvilinear fibre composite laminated panels, J. Compos. Mater. 48 (22) (2014) 2761–2777,
[Link]
[5] C.A. Shankara, N.G.R. Iyengar, A C0 element for the free vibration analysis of laminated
Pr

composite plates, J. Sound Vib. 191 (5) (1996) 721-738,

[Link]

41

This preprint research paper has not been peer reviewed. Electronic copy available at: [Link]
 [6] M. Boscolo, J.R. Banerjee, Layer-wise dynamic stiffness solution for free vibration analysis of
laminated composite plates, J. Sound Vib. 333 (1) (2014) 200–227.

ed
[Link]
[7] H. Akhavan, P. Ribeiro, Natural modes of vibration of variable stiffness composite laminates
with curvilinear fibers, Compos. Struct. 93 (11) (2011) 3040–3047,
[Link]

iew
[8] S. Honda, Y. Narita, Natural frequencies and vibration modes of laminated composite plates
reinforced with arbitrary curvilinear fiber shape paths, J. Sound Vib. 331 (1) (2012) 180–191,
[Link]
[9] Y.H. Lim, Finite-element simulation of closed loop vibration control of a smart plate under
transient loading, Smart Materials and Structures,12 (2) (2003) 272,
[Link]

v
[10] S.B. Kerur, A. Ghosh, Active vibration control of composite plate using AFC actuator and
PVDF sensor, International Journal of Structural Stability and Dynamics, 11 (2) (2011) 237–

re
255, [Link]
[11] N. Sharma, P.K. Swain, D.K. Maiti, B.N. Singh, Static and free vibration analyses and
dynamic control of smart variable stiffness laminated composite plate with delamination,
Compos. Struct. 280 (2022) 114793, [Link]
er
[12] P.K. Mahato, D.K. Maiti, Flutter Control of Smart Composite Structures in Hygrothermal
Environment, Journal of Aerospace Engineering 23 (4) (2010) 317-326,
[Link]
pe
[13] J. Duan, B. Xu, X. Xue, L. Shi, P. Hao, Subsonic aeroelastic behaviors of a composite plate
with embedded MFC actuators under hygrothermal environment, Thin-Walled Structures 200
(2024) 111906, [Link]
[14] P. Mondal, J.P. Varun, P.K. Mahato, Open loop flutter control of optimally oriented smart
variable stiffness plates under hygrothermal environment, European Journal of Mechanics,
A/Solids 106 (2024) 105284, [Link]
ot

[15] N. Sharma, P.K. Swain, D.K. Maiti, Active flutter suppression of damaged variable stiffness
laminated composite rectangular plate with piezoelectric patches, Mechanics of Advanced
Materials and Structures 31 (6) (2024) 1229–1249,
tn

[Link]
[16] S. Raja, A.A. Pashilkar, R. Sreedeep, J. V. Kamesh, Flutter control of a composite plate with
piezoelectric multilayered actuators, Aerosp. Sci. Technol. 10 (5) (2006) 435–441,
[Link]
rin

[17] Z.G. Song, F.M. Li, E. Carrera, P. Hagedorn, A new method of smart and optimal flutter
control for composite laminated panels in supersonic airflow under thermal effects, J. Sound
Vib. 414 (2018) 218–232, [Link]
[18] J.A. Moreira, F. Moleiro, A.L. Araújo, A. Pagani, Active aeroelastic flutter control of
ep

supersonic smart variable stiffness composite panels using layerwise models, Compos. Struct.
343 (2024) 118287, [Link]
[19] J.A. Moreira, F. Moleiro, A.L. Araújo, A. Pagani, Active aero-visco-elastic flutter control and
layerwise modelling of supersonic smart sandwich panels with variable stiffness composites,
Pr

Aerosp. Sci. Technol. 157 (2025) 109847, [Link]

[20] R.M. Orris, M. Petyt, A finite element study of the vibration of trapezoidal plates, J. Sound
Vib. 27 (3) (1973) 325–344, [Link]
42

This preprint research paper has not been peer reviewed. Electronic copy available at: [Link]
 [21] R.S. Srinivasan, B.J.C. Babu, Free vibration of cantilever quadrilateral plates, The Journal of
the Acoustical Society of America, 73 (3) (1983) 851-855, [Link]

ed
[22] R.S. Srinivasan, B.J.C. Babu, Flutter analysis of cantilevered quadrilateral plates, J. Sound
Vib. 98(1) (1985) 45-53, [Link]
[23] R.S. Srinivasan, B.J.C. Babu, Free vibration and flutter of laminated quadrilateral plates,
Computers & structures, 27 (2) (1987) 297-304, [Link]

iew
6.
[24] H.T.Y. Yang, Y.C. Wu, A geometrically non‐linear tensorial formulation of a skewed
quadrilateral thin shell finite element, Int. J. Numer. Methods Eng. 28 (12) (1989) 2855–2875,
[Link]
[25] R.M.V Pidaparti, C.C. Chang, Finite element supersonic futter analysis of skewed and cracked

v
composite panels, Computers & structures 69 (2) (1998) 265-270,
[Link]

re
[26] T.V.R. Chowdary, P.K. Sinha, S. Parthan, Finite element flutter analysis of composite skew
panels, Computers & structures 58 (3) (1996) 613-620, [Link]
7949(95)00153-8.
[27] M. Barik, M. Mukhopadhyay, Finite element free flexural vibration analysis of arbitrary
er
plates, Finite elements in analysis and design 29 (2) (1998) 137-151,
[Link]
[28] M.K. Singha, M. Ganapathi, Large amplitude free flexural vibrations of laminated composite
pe
skew plates, Int. J. Non Linear Mech. 39 (10) (2004) 1709–1720,
[Link]
[29] M.K. Singha, M. Ganapathi, A parametric study on supersonic flutter behavior of laminated
composite skew flat panels, Compos. Struct. 69 (1) (2005) 55–63,
[Link]
ot

[30] R.M. Kanasogi, M.C. Ray, Control of geometrically nonlinear vibrations of skew laminated
composite plates using skew or rectangular 1-3 piezoelectric patches, International Journal of
Mechanics and Materials in Design 9 (2013) 325–354, [Link]
tn

9224-z.
[31] K.M. Liew, Vibration of symmetrically laminated cantilever trapezoidal composite plates, Int.
J. Mech. Sc. 34 (4) (1992) 299-308. [Link]
rin

[32] S. Wang, Vibration of thin skew fibre reinforced composite laminates, J. Sound Vib. 201 (3)
(1997) 335-352, [Link]
[33] A. Houmat, Nonlinear free vibration analysis of variable stiffness symmetric skew laminates,
European Journal of Mechanics, A/Solids 50 (2015) 70–75,
ep

[Link]
[34] V. Khalafi, J. Fazilati, Supersonic panel flutter of variable stiffness composite laminated skew
panels subjected to yawed flow by using NURBS-based isogeometric approach, J. Fluids
Struct. 82 (2018) 198–214, [Link]
Pr

[35] T. Farsadi, D. Asadi, H. Kurtaran, Fundamental frequency optimization of variable stiffness

composite skew plates, Acta Mech. 232 (2021) 555–573, [Link]
02871-9.

43

This preprint research paper has not been peer reviewed. Electronic copy available at: [Link]
 [36] M. Rahmanian, T. Farsadi, H. Kurtaran, Nonlinear flutter of tapered and skewed cantilevered
plates with curvilinear fiber paths, J. Sound Vib. 500 (2021) 116021,

ed
[Link]
[37] K. Torabi, H. Afshari, F.H. Aboutalebi, Vibration and flutter analyses of cantilever trapezoidal
honeycomb sandwich plates, Journal of Sandwich Structures and Materials 21 (8) (2019)
2887–2920, [Link]

iew
[38] P.K. Mahato, Static, Dynamic and Flutter Control of Laminated Composite Plates in
Hygrothermal Environment Employing Active Fiber Composite,PhD Thesis, IIT Kharagpur,
(2010), [Link]
[39] S. Raja, P.K. Sinha, G. Prathap, Active Stiffening and Active Damping Effects on Closed
Loop Vibration Control of Composite Beams and Plates, Journal of reinforced plastics and

v
composites 22 (12) (2003) 1101-1121, [Link]
[40] K.J. Bathe, Finite Element Procedures, 2nd ed., Prentice Hall, (1996).

re
[41] Aeroelastic Analysis User’s Guide, MSC Software Corporation, USA, (2019).

er
pe
ot
tn
rin
ep
Pr

44

This preprint research paper has not been peer reviewed. Electronic copy available at: [Link]