Aeroelasticity and dynamic instability research of carbon nanotube reinforced composite laminated cylindrical shells under hygrothermal environment

Abstract

This study presents the aerodynamic stability and dynamic instability of hierarchical cylindrical shells reinforced with carbon nanotubes (CNTs) and surface-bonded piezoelectric composite films. The analysis employs double-ended free boundary conditions under hygrothermal environment, explicitly accounting for the relationship between material properties and hygrothermal effects. Fundamental kinetic equations of the CNT-reinforced, piezoelectrically layered cylindrical shells are derived from Love’s classical thin-shell theory, which excludes the pyroelectric effects of the piezoelectric layers. The Rayleigh–Ritz method is employed for the numerical solution of the dynamic equations, through combining with Bolotin’s approach. Chebyshev polynomials discretize the displacement fields. In the meanwhile, aerodynamic loading is characterized using first-order linear piston theory with airflow deflection. Furthermore, the analysis systematically examines the influence of critical parameters, including Mach number, airflow deflection angle, CNT volume fraction, fiber orientation angle, temperature, humidity, static load parameters, piezoelectric layer properties, and applied voltage, on the aerodynamic stability boundaries and regions of dynamic instability.

Keywords

free-free cylindrical shells aeroelasticity dynamic instability piezoelectric carbon nanotubes

Introduction

Composite-layered cylindrical shell has been widely used in the field of diverse specializations, ranging from mechanical, civil, aerospace, to architectural engineering. The composite hierarchical cylindrical shells reinforced with CNTs has attracted considerable interest in aerospace applications due to their superior material properties and design flexibility. These structures are frequently exposed to complex loading conditions, such as aerodynamic, thermal, and hygroscopic effects. Consequently, a critical investigation of CNT-reinforced hierarchical cylindrical shells is essential in aerospace applications.

In recent years, CNTs have been increasingly employed to reinforce the matrix phase of composite materials, owing to their excellent specific strength and stiffness. However, CNT composite materials experience significant sensitivity to hygrothermal effects, especially under combined hygrothermal-mechanical loading. As shown by Adams and Miller (1977), Bowles and Tompkins (1989), and Boukhoulda et al. (2006), such environmental exposure progressively degrades their elastic modulus and overall mechanical strengths.

Contemporary studies demonstrate the problem through the lens of CNT distribution, which are typically assuming either homogeneous dispersion or functionally graded (FG) patterning within the matrix, and idealized boundary conditions, such as commonly simply-supported (S-S) or clamped-clamped (C-C) at the ends. For instance, the work by Arani et al. (2019) used Third-Order Shear Deformation Theory (TSDT) coupled with linear piston theory to investigate how FG-CNT material gradation, yaw angle, and geometry dictate aerodynamic stability boundaries. This naturally leads to considerations of boundary constraint effects, as compared for aeroelastic stability by Asadi and Wang (2017) through S-S and C-C configurations. Further probing the dynamic consequences of supersonic flow, Avramov et al. (2019) examined the nonlinear self-excited vibration in S-S FG-CNT shells with the harmonic balance method. Chakraborty et al. (2021) further analyzed the parametric instability of open FG-CNT shells under periodic loading, incorporating Rayleigh structural damping. While, Heydarpour and Malekzadeh (2018) quantified the influence of rotational speed and FG-CNT gradation on the instability regions. In further, Mohammadi et al. (2021) tackled the challenging case of axially moving shells with time-varying velocity using multiscale analysis. Efforts on control are also made by Zhang et al. (2017) who explored nonconstant aerodynamic stability and demonstrated flutter suppression strategies in S-S FG-CNT shells using the principle of energy diversity. In terms of the thermal buckling of cylindrical shells, Heidari-Soureshjani et al. (2024a) established a sophisticated multi-field analytical framework to reveal novel buckling vibration characteristics of porous FG conical shells under hygrothermal coupling, and explicitly quantify the coupled degradation effects of temperature, humidity, and porosity on both vibration frequency and buckling strength. To expand the structural and thermal analysis, Heidari-Soureshjani et al. (2024b) subsequently investigated the dynamic stability of rotating, high-temperature FG-CNT biconical shells, and employed FSDT and the GDQ method for studying the roles of temperature, rotation speed, joining parameters, and CNT properties. The findings on complex geometries under thermal loads are done by Pakravan et al. (2022), who provided a characterization of the fundamental free vibration and buckling mechanisms governing CNT-reinforced conical shells in thermal environments, based on Haar wavelets and FSDT. Besides, the work also details the influence of geometrical dimensions, boundary conditions, temperatures, and material properties.

While significant research exists on the flutter problem in CNT-reinforced functionally graded cylindrical shells under simply supported and clamped boundary conditions, the flutter and dynamic instability boundaries for multi-layered CNT-reinforced cylindrical shells under free-free (F-F) boundary conditions remains unaddressed. This gap is critical given the extensive use of F-F cylindrical shells in aerospace structures, , such as the intercontinental ballistic missile bodies and rocket engine thermal insulation liners, where high-speed flight subjects them to supersonic airflow and thrust oscillations (Li et al., 2018). Addressing this need, the present study introduces two key advancements. Firstly, the effective material properties of the CNT-reinforced multi-layered cylindrical shells are characterized using the Mori-Tanaka micromechanical method, departing from the less accurate mixture rule commonly employed in prior FG-CNT studies. Secondly, we derive discrete flutter and dynamic instability models for F-F boundary conditions by directly applying the Ritz minimization procedure to the Lagrange energy functional. The proposed methodology and resulting findings have been rigorously validated against established literature. Furthermore, aeroelastic stability boundaries and dynamic instability for F-F hierarchical shells under supersonic flow conditions are systematically studied.

Numerical approach and model development

Figure 1(a) illustrates the geometrical size of a CNT-reinforced composite cylindrical shell structure with surface-bonded piezoelectric actuator layers on both inner and outer surfaces. Boundary conditions are implemented via an artificial spring methodology, which represent schematically at the shell edges with associated stiffness coefficients. The $x - θ$ plane, marked by the dashed line, signifies the midplane equally distant from top and bottom surfaces of the shells. $L$ represents the length of the shells, while $R$ represents distance from the center of the circle to the mid-plane. The overall thickness of the primary composite shells, $h$ , is calculated as the sum of the thicknesses of each layer. Furthermore, $h_{p}$ represents the thickness of the piezoelectric plies. The coordinate system defines $x$ as the axial direction, $θ$ as the circumferential direction, and $z$ as the normal direction relative to the midplane. As depicted in Figure 1(b), axes $1, 2, 3$ represent the material’s principal longitudinal, transverse, and normal directions, respectively, where $Θ$ denotes the counterclockwise angle between the shell’s axial line ( $x -$ axis) and the CNT fiber axis ( $1 -$ axis).

Figure 1.

(a) CNT reinforced multilayer cylindrical shells with piezoelectric layers in supersonic flow subject to spring constraints and local enlargement diagram. (b) Single-layer linear CNT reinforced composite shells and coordinates.

Material parameter

We employ the Mori–Tanaka micromechanical model to study the effective material properties of CNT-reinforced cylindrical shells. Specifically, straight CNTs are embedded in the cylindrical shell matrix to enhance mechanical performance, with the horizontal orientation angle $Θ$ defined between the CNT fibers and the shell’s symmetry axis. The Hill’s elastic modulus of the composite can be expressed as (Hill, 1965)

k = \frac{E_{m} {E_{m} + f_{m} + 2 k_{r} (1 + v_{m}) [1 + f_{r} (1 - 2 v_{m})]}}{2 (1 + v_{m}) [E_{m} (1 + f_{r} - 2 v_{m}) + 2 f_{m} k_{r} (1 - v_{m} - 2 v_{m}^{2})]}

(1)

l = \frac{E_{m} {f_{m} v_{m} [E_{m} + 2 k_{r} (1 + v_{m})] + 2 f_{r} k_{r} (1 - v_{m}^{2})}}{(1 + v_{m}) [2 f_{m} k_{r} (1 - v_{m} - 2 v_{m}^{2}) + E_{m} (1 + f_{r} - 2 v_{m})]}

(2)

\begin{array}{l} n = & \frac{E_{m}^{2} f_{m} (1 + f_{r} - f_{m} v_{m}) + 2 f_{r} f_{m} (k_{r} n_{r} - l_{r}^{2}) {(1 + v_{m})}^{2} (1 - 2 v_{m})}{(1 + v_{m}) [2 f_{m} k_{r} (1 - v_{m} - 2 v_{m}^{2}) + E_{m} (1 + f_{r} - 2 v_{m})]} \\ + \frac{E_{m}^{2} [2 f_{m}^{2} k_{r} (1 - v_{m}) + 2 n_{r} f_{m} (1 - 2 v_{m} + f_{r}) + 4 f_{m} f_{r} l_{r} v_{m}]}{2 f_{m} k_{r} (1 - v_{m} - 2 v_{m}^{2}) + E_{m} (1 + f_{r} - 2 v_{m})} \end{array}

(3)

p = \frac{E_{m} [E_{m} f_{m} + 2 (1 + f_{r}) p_{r} (1 + v_{m})]}{2 (1 + v_{m}) [E_{m} (1 + f_{r}) + 2 f_{m} p_{r} (1 + v_{m})]}

(4)

m = \frac{E_{m} [E_{m} f_{m} + 2 m_{r} (1 + v_{m}) (3 + f_{r} - 4 v_{m})]}{2 (1 + v_{m}) {E_{m} [f_{m} + f_{r} (1 - v_{m})] + 2 f_{m} m_{r} (3 - v_{m} - 4 v_{m}^{2})}}

(5)

where in above equations,

k_{r}

l_{r}

n_{r}

p_{r}

m_{r}

represent the Hill’s modulus of elasticity of CNT fibers,

E_{m}

and

v_{m}

represent the modulus of elasticity and Poisson’s ratio of the matrix,

f_{r} = 1 - f_{m}

means the percentage of carbon nanotubes in the matrix, that is, the volume fraction. Accordingly, the effective parameters of the material of CNT reinforced composite, including elastic modulus, modulus of shear and Poisson’s ratio, can be given as (Bisheh et al., 2020)

E_{11} = b - \frac{l^{2}}{k}, E_{22} = \frac{4 m (k n - l^{2})}{k n - l^{2} + m n}, G_{12} = p, v_{12} = \frac{1}{2 k}

(6)

Hydrothermal expansion coefficient of composite’s plies

According to Bisheh et al. (2020), the principal longitudinal ( $α_{11}$ ) and transverse ( $α_{22}$ ) thermal expansion coefficients (CTEs) for composite material are described as

α_{11} = \frac{E_{11 r} α_{11 r} f_{r} + E_{m} α_{m} f_{m}}{E_{11 r} f_{r} + E_{m} f_{m}}, α_{22} = α_{22 r} f_{r} (1 + v_{12 r} \frac{α_{11 r}}{α_{22 r}}) + (1 + v_{m}) α_{m} f_{m} - (v_{12 r} f_{r} + v_{m} f_{m}) α_{11}

(7)

where the elastic modulus of CNTs and matrix can be written as

E_{11 r}

and

E_{m}

. The coefficients of thermal expansion (CTEs) for the matrix phase is indicated by

α_{m}

. For the reinforcing fibers, the longitudinal and transverse orientations CTEs are represented by

α_{11 r}

and

α_{22 r}

, respectively. Additionally,

v_{m}

and

v_{12 r}

denote the Poisson’s ratios for the matrix and the reinforcing fibers, correspondingly. Ignoring coefficients of moisture expansion (CMEs) of fibers (Bisheh et al., 2020), the principle longitudinal and transverse CMEs for the fiber-reinforced composite are given as

β_{11} = β_{m} \frac{E_{m} f_{m}}{E_{11 r} f_{r} + E_{m} f_{m}}, β_{22} = β_{m} \frac{f_{m}}{E_{11 r} f_{r} + E_{m} f_{m}} [(E_{11 r} f_{r} + E_{m} f_{m}) + f_{r} (v_{m} E_{11 r} - v_{12 r} E_{m})]

(8)

where

β_{m}

is the CMEs of the matrix.

To assess the thermomechanical properties of CNTs as a function of temperature, a third-order interpolation method is employed, as detailed below (Dat et al., 2020)

E_{11 r} (T) [T P a] = 6.3998 - 4.33817 \times 1 0^{- 3} T + 7.34 \times 1 0^{- 6} T^{2} - 4.45833 \times 1 0^{- 9} T^{3}

(9)

α_{11 r} (T) [\frac{1 0^{- 6}}{K}] = - 1.12515 + 0.02291688 T - 2.887 \times 1 0^{- 5} T^{2} + 1.13625 \times 1 0^{- 8} T^{3}

(10)

α_{22 r} (T) [\frac{1 0^{- 6}}{K}] = 5.43715 - 0.984625 T + 2.9 \times 1 0^{- 7} T^{2} - 1.25 \times 1 0^{- 11} T^{3}

(11)

E_{m} (T) [G P a] = (3.52 - 0.0034 Δ T) G P a

(12)

α_{m} (T) [\frac{1 0^{- 6}}{K}] = 45 (1 + 0.0005 Δ T), β_{m} = 0.06

(13)

In calculation of equations (9)–(11)

T = T_{0} + Δ T

, where

T_{0} = 300 K

is the reference temperature. The CTEs and CMEs of composite structures are described by

α_{11}

α_{22}

β_{11}

and

β_{22}

, these coefficients can be transformed into general coordinates as follow

α_{x} = α_{11} \cos^{2} (Θ) + α_{11} \sin^{2} (Θ), α_{θ} = α_{11} \sin^{2} (Θ) + α_{22} \cos^{2} (Θ), α_{x θ} = 2 (α_{11} - α_{22}) \sin (Θ) \cos (Θ)

(14)

β_{x} = β_{11} \cos^{2} (Θ) + β_{11} \sin^{2} (Θ), β_{θ} = β_{11} \sin^{2} (Θ) + β_{22} \cos^{2} (Θ), β_{x θ} = 2 (β_{11} - β_{22}) \sin (Θ) \cos (Θ)

(15)

Constitutive relation and energy

The Love classical shells theory indicates the relationship between its mid-surface displacement and strain as follows

ε_{0} = {[\frac{\partial u}{\partial x}, \frac{1}{R} (\frac{\partial u}{\partial θ} + w), \frac{\partial v}{\partial x} + \frac{1}{R} \frac{\partial u}{\partial θ}]}^{T}

(16)

κ = {[- \frac{\partial^{2} w}{\partial x^{2}}, - \frac{1}{R^{2}} (\frac{\partial^{2} w}{\partial θ^{2}} - \frac{\partial v}{\partial θ}), - \frac{2}{R} (\frac{\partial^{2} w}{\partial x \partial θ} - \frac{\partial v}{\partial x})]}^{T}

(17)

where

ε_{0} = {[ε_{x}^{0}, ε_{θ}^{0}, ε_{x θ}^{0}]}^{T}

κ = {[χ_{x}, χ_{θ}, χ_{x θ}]}^{T}

are components of membrane strain and curvature variation of the curved surface in the shells with cylindrical geometry, correspondingly. The components of the displacement of the coordinate system at arbitrary points in the mid-plane can be characterized as

u

v

, and

w

. The components of strain field at arbitrary points on the shells can be summarized as

ε = ε_{0} + z κ

(18)

where

ε = [ε_{x}, ε_{θ}, ε_{x θ}]

represents the sector of mechanical in-plane deformations at any location inside the cylindrical shells, which depends on the reference surface in-plane strains, the vertical distance

z

from this surface, and the shells’ in-plane curvatures.

When subjected to mechanical and hygrothermal loads, the in-plane stress distributions within the laminate can be mathematically described as follows

{\begin{cases} σ_{x} \\ σ_{θ} \\ σ_{x θ} \end{cases}}_{K}^{M} = T_{K}^{- 1} [\begin{array}{l} Q_{11} & Q_{12} \\ Q_{12} & Q_{22} \\ Q_{66} \end{array}] {[T_{K}^{- 1}]}^{T} {\begin{cases} ε_{x} \\ ε_{θ} \\ ε_{x θ} \end{cases}}, {\begin{cases} σ_{x} \\ σ_{θ} \\ σ_{x θ} \end{cases}}_{K}^{T} = - T_{K}^{- 1} [\begin{array}{l} Q_{11} & Q_{12} \\ Q_{12} & Q_{22} \\ Q_{66} \end{array}] {[T_{K}^{- 1}]}^{T} {\begin{cases} α_{x} \\ α_{θ} \\ α_{x θ} \end{cases}} Δ T, {\begin{cases} σ_{x} \\ σ_{θ} \\ σ_{x θ} \end{cases}}_{K}^{C} = - T_{K}^{- 1} [\begin{array}{l} Q_{11} & Q_{12} \\ Q_{12} & Q_{22} \\ Q_{66} \end{array}] {[T_{K}^{- 1}]}^{T} {\begin{cases} β_{x} \\ β_{θ} \\ β_{x θ} \end{cases}} Δ C

(19)

Q_{11} = \frac{E_{11}}{1 - ν_{12} ν_{21}}, Q_{12} = \frac{ν_{21} E_{11}}{1 - ν_{12} ν_{21}}, Q_{21} = \frac{ν_{21} E_{11}}{1 - ν_{12} ν_{21}}, Q_{22} = \frac{E_{22}}{1 - ν_{12} ν_{21}}, Q_{66} = G_{12}

(20)

where superscript

M

T

, and

C

stand for mechanical, thermal, and moisture stress, correspondingly.

Δ T

and

Δ C

represent the variations in temperature and humidity. The coordinate transformation matrix

T_{k}

as follows

T_{k} = [\begin{array}{c} \cos^{2} (Θ) & \sin^{2} (Θ) & \sin (Θ) \cos (Θ) \\ \sin^{2} (Θ) & \cos^{2} (Θ) & - \sin (Θ) \cos (Θ) \\ - 2 \sin (Θ) \cos (Θ) & 2 \sin (Θ) \cos (Θ) & \cos^{2} (Θ) - \sin^{2} (Θ) \end{array}]

(21)

The three stress vectors in equation (19) can be, respectively, abbreviated as

σ_{K}^{M} = {\bar{Q}}_{K} ε_{K}, σ_{K}^{M} = - {\bar{Q}}_{K} α_{K} Δ T, σ_{K}^{M} = - {\bar{Q}}_{K} β_{K} Δ C

(22)

By summing the integrals of stress over each thickness of the shells can evaluate the internal force

N = {[N_{x} N_{θ} N_{x θ}]}^{T}

and moment

M = {[M_{x} M_{θ} M_{x θ}]}^{T}

N = \sum_{K = 1}^{N} \int_{z_{k - 1}}^{z_{k}} σ_{K}^{M} d z, M = \sum_{K = 1}^{N} \int_{z_{k - 1}}^{z_{k}} z σ_{K}^{M} d z

(23)

where

N

is used to represents the aggregate number of composite plies, whereas

z_{k - 1}

and

z_{k}

denote the

z -

coordinates of the

K^{t h}

layer, respectively. By substituting equation (23) into equation (22) and integrating over the aggregate thickness of the main shells yields the following results

N = A ε_{0} + B κ

(24)

M = B ε_{0} + D κ

(25)

where A,

B

and

D

are the stiffness matrix of membrane, stiffness matrix of coupling and stiffness matrix of bending of the composite cylindrical, correspondingly, and

A = \sum_{K - 1}^{N} {\bar{Q}}_{K} (z_{K} - z_{K - 1}), B = \frac{1}{2} \sum_{K - 1}^{N} {\bar{Q}}_{K} (z_{K}^{2} - z_{K - 1}^{2}), D_{K} = \frac{1}{3} \sum_{K - 1}^{N} {\bar{Q}}_{K} (z_{K}^{3} - z_{K - 1}^{3}),

(26)

Similarly,

N_{Δ T}

and

M_{Δ T}

are used to represent the force and moment resultant vectors induced by thermal effects in the host substrate laminated shells, correspondingly.

N_{Δ C}

M_{Δ C}

are the force and moment tensor caused by humidity-induced of the host substrate hierarchical shells, and they can be characterized as

N_{Δ T} = \sum_{K - 1}^{N} {\bar{Q}}_{K} (z_{K} - z_{K - 1}) α_{K} Δ T, M_{Δ T} = \sum_{K - 1}^{N} {\bar{Q}}_{K} (z_{K}^{2} - z_{K - 1}^{2}) α_{K} Δ T

(27)

N_{Δ C} = \sum_{K - 1}^{N} {\bar{Q}}_{K} (z_{K} - z_{K - 1}) β_{K} Δ C, M_{Δ C} = \sum_{K - 1}^{N} {\bar{Q}}_{K} (z_{K}^{2} - z_{K - 1}^{2}) β_{K} Δ C

(28)

where

N_{Δ T} = {[N_{T x} N_{T θ} N_{T x θ}]}^{T}

N_{Δ C} = {[N_{C x} N_{C θ} N_{C x θ}]}^{T}

The energy of the hierarchical shells structure, which encompasses the strain energy as well as the energy of spring and kinetic, can be characterized as follows

U = \frac{1}{2} \iint_{S} [N^{T} ε_{0} + M^{T} κ] d S + \frac{1}{2} \int_{0}^{2 π} [k_{u} u^{2} + k_{v} v^{2} + k_{w} w^{2} + \frac{K_{w}}{L^{2}} {(\frac{\partial w}{\partial x})}^{2}] R d θ

(29)

T_{c} = \frac{1}{2} \iint_{S} ρ_{T} ({\dot{u}}^{2} + {\dot{v}}^{2} + {\dot{w}}^{2}) d S

(30)

where

ρ_{T}

is computed by integrating

ρ_{r}

(carbon nanotubes) and

ρ_{m}

(matrix)

ρ_{T} = \int_{- h / 2}^{h / 2} (ρ_{r} f_{r} + ρ_{m} f_{m}) d z

(31)

The work done by the hygrothermal loads are expressed as (Wang et al., 2023)

W_{Δ T} = - \frac{1}{2} \iint_{S} [N_{T x} {(\frac{\partial w}{\partial x})}^{2} + N_{T θ} {(\frac{\partial w}{R \partial θ})}^{2} + 2 N_{T x θ} \frac{\partial w}{\partial x} \frac{\partial w}{R \partial θ}] R d S

(32)

W_{Δ C} = - \frac{1}{2} \iint_{S} [N_{C x} {(\frac{\partial w}{\partial x})}^{2} + N_{C θ} {(\frac{\partial w}{R \partial θ})}^{2} + 2 N_{C x θ} \frac{\partial w}{\partial x} \frac{\partial w}{R \partial θ}] R d S

(33)

The energy due to the application of an axial compressive force can be represented as (Asadi and Wang, 2017)

U_{g} = \iint_{S} N_{x} [\frac{\partial u}{\partial x} - \frac{1}{2} {(\frac{\partial^{2} w}{\partial x^{2}})}^{2}] R d S

(34)

where

N_{x}

is the axial compressive load. The work performed by the aerodynamic load can be defined as (Hasheminejad and Aghayi Motaaleghi, 2015)

W_{Δ P} = - \iint_{S} \frac{γ p_{\infty} M_{a}^{2}}{\sqrt{M_{a} - 1}} (\frac{\partial w}{\partial x} \cos Λ + \frac{1}{R} \frac{\partial w}{\partial θ} \sin Λ + \frac{1}{M_{a} a_{\infty}} (\frac{M_{a}^{2} - 2}{M_{a}^{2} - 1}) \dot{w} - \frac{1}{2 R \sqrt{M_{a} - 1}} w) w R d S

(35)

where

γ

p_{\infty}

a_{\infty}

M_{a}

Λ

are the air adiabatic index, the static pressure due to the air flow, the local sonic velocity, the Mach number and the airflow deflection angle, respectively.

Constitutive relation of piezoelectric film

The piezoelectric film, having a thickness denoted as $h_{p}$ , are adhered to the external and internal edges of the primary hierarchical shells, as shown in Figure 1(b). Polarization is achieved by applying an external electric field in the thickness direction. The coupled electromechanical constitutive equation for a piezocomposite layer is presented as follows (Chai et al., 2017)

σ^{p} = {\bar{c}}^{p} ε^{p} - \bar{e} E

(36)

D = {\bar{e}}^{T} ε^{p} + Ξ E

(37)

The above condensed form of equations further expanded as below

{\begin{cases} σ_{x x} \\ σ_{θ θ} \\ σ_{x θ} \end{cases}}^{p} = {[\begin{array}{l} {\bar{c}}_{11} & {\bar{c}}_{12} & {\bar{c}}_{16} \\ {\bar{c}}_{22} & {\bar{c}}_{26} \\ s y s & {\bar{c}}_{66} \end{array}]}^{p} {\begin{cases} ε_{x x} \\ ε_{θ θ} \\ ε_{x θ} \end{cases}}^{p} - [\begin{array}{l} 0 & 0 & {\bar{e}}_{31} \\ 0 & 0 & {\bar{e}}_{32} \\ 0 & 0 & {\bar{e}}_{33} \end{array}] {\begin{cases} 0 \\ 0 \\ E_{z} \end{cases}}

(38)

{\begin{cases} 0 \\ 0 \\ D_{z} \end{cases}} = [\begin{array}{l} 0 & 0 & {\bar{e}}_{31} \\ 0 & 0 & {\bar{e}}_{32} \\ 0 & 0 & {\bar{e}}_{33} \end{array}] {\begin{cases} ε_{x x} \\ ε_{θ θ} \\ ε_{x θ} \end{cases}}^{p} - [\begin{array}{l} 0 & 0 & 0 \\ 0 & 0 & 0 \\ 0 & 0 & {\bar{Ξ}}_{33} \end{array}] {\begin{cases} 0 \\ 0 \\ E_{z} \end{cases}}

(39)

where the superscript

p

indicates that the stress and strain parameters in above equations correspond to the physical properties of the piezoelectric fibers. The tensor

{σ}^{p}

represents the stress, the tensor

{ε}^{p}

represents the strain, the stiffness matrix, once transformed and reduced, is identified as

{[\bar{c}]}^{p}

, the matrix

[\bar{e}]

represents the piezoelectric constant, the tensor

[\bar{Ξ}]

represents the dielectric constant, the tensor

[E]

depend with voltage and the tensor

[D]

denotes the electric displacement of an piezoelectric material. The angle-ply of the piezoelectric fibers is defined as, the transformation matrix is also based on equation (21), but with the variable

Θ

replaced by the variable.

{[\bar{c}]}^{p}

and

[\bar{e}]

can be expressed as

{[\bar{c}]}^{p} = T_{K}^{- 1} [\begin{array}{l} c_{11}^{p} & c_{12}^{p} \\ c_{12}^{p} & c_{22}^{p} \\ c_{66}^{p} \end{array}] {[T_{K}^{- 1}]}^{T}, [\bar{e}] = T_{K}^{- 1} [\begin{array}{l} 0 & 0 & 0 \\ 0 & 0 & 0 \\ e_{31} & e_{32} & 0 \end{array}]

(40)

c_{11}^{p} = \frac{E_{11}^{p}}{1 - ν_{12}^{p} ν_{21}^{p}}, c_{12}^{p} = \frac{ν_{21}^{p} E_{11}^{p}}{1 - ν_{12}^{p} ν_{21}^{p}}, c_{21}^{p} = \frac{ν_{12}^{p} E_{12}^{p}}{1 - ν_{12}^{p} ν_{21}^{p}}, c_{22}^{p} = \frac{E_{22}^{p}}{1 - ν_{12}^{p} ν_{21}^{p}}, c_{66}^{p} = G_{12}^{p}

(41)

where

E_{11}^{p}

E_{22}^{p}

G_{12}^{p}

ν_{12}^{p}

and

ν_{21}^{p}

are be used to mean the modulus of elastic, modulus of shear, and Poisson’s ratio associated with the piezoelectric fibers, correspondingly. Further, the force resultant vector

N^{p}

, moment resultant vector

M^{p}

of the piezoelectric layer can be represented as

N^{p} = A^{p} ε_{0} + B^{p} κ

(42)

M^{p} = B^{p} ε_{0} + D^{p} κ

(43)

where

A^{p}

B^{p}

, and

D^{p}

are the matrix of membrane stiffness, matrix of coupling stiffness and matrix of bending stiffness of the piezoelectric film, correspondingly.

[A^{p}, B^{p}, D^{p}] = \int_{h / 2}^{h / 2 + h_{p}} {[\bar{c}]}^{p} [1, z, z^{2}] d z + \int_{- h / 2 - h}^{h / 2_{p}} {[\bar{c}]}^{p} [1, z, z^{2}] d z

(44)

Only the outer piezoelectric film is used as the actuator,

N_{V}^{p}

and

M_{V}^{p}

represent the force and moment resultant vectors, respectively, induced by voltage in the outer piezoelectric layer

[N_{V}^{p}, M_{V}^{p}] = \int_{h / 2}^{h / 2 + h_{p}} - \bar{e} E [1, z] d z

(45)

The energy of strain and energy of kinetic associated with the inner and outer piezoelectric film can be formulated as

U_{p} = \frac{1}{2} \iint_{S} [N^{p, T} ε_{0} + M^{p, T} κ] d S

(46)

T_{p} = \frac{1}{2} \iint_{S} ρ_{p} ({\dot{u}}^{2} + {\dot{v}}^{2} + {\dot{w}}^{2}) d S

(47)

where

ρ_{p}

represents the mass surface density of the piezoelectric film. The energy done by the voltage applied on the outer surface piezoelectric film is attained as (Zhou et al., 2019)

W_{E} = - \frac{1}{2} \iint_{S} [N_{V, x}^{p} {(\frac{\partial w}{\partial x})}^{2} + N_{V, θ}^{p} {(\frac{\partial w}{\partial θ})}^{2} + 2 N_{V, x θ}^{p} (\frac{\partial w}{\partial x} \frac{\partial w}{\partial θ})] R d S

(48)

Governing equations

The displacement field of a multi-layered cylindrical shell can be mathematically expressed in terms of axial orthogonal polynomials combined with circumferential Fourier series

u (η, θ, t) = \sum_{m = 1}^{M} U_{m} P_{m} \cos (n θ) e^{i Ω t}, v (η, θ, t) = \sum_{m = 1}^{M} V_{m} P_{m} \sin (n θ) e^{i Ω t}, w (η, θ, t) = \sum_{m = 1}^{M} W_{m} P_{m} \cos (n θ) e^{i Ω t}

(49)

where the unknown coefficients can be named by

U_{m}

V_{m}

, and

W_{m}

, the angular frequency can be denoted by

Ω

,the mode number of the shells in the

x

and

θ

directions can be expressed by

m

and

n

, the regular expression for the first kind of Chebyshev polynomials, denoted as

p_{m} (η)

, can be described as follow

P_{1} (η) = \frac{a + b + 2}{2} φ - \frac{a - b}{2}, {P_{0} (η) = 1, P}_{m} (η) = \frac{(a + b + 2 m - 1) [a^{2} - b^{2} + φ (a + b + 2 m) (a + b + 2 m - 2)]}{2 m (a + b + m) (a + b + 2 m - 2)} {P_{m} (η - 1)}_{a, b} = - \frac{(a + m - 1) (b + m - 1) (a + b + 2 m)}{2 m (a + b + m) (a + b + 2 m - 2)} P_{m - 2} (η - 1) = - 1 / 2; m = 2, 3 \cdot \cdot \cdot

(50)

Since the Chebyshev polynomial is defined within the range of

[- 1, 1]

, this study necessitates a transformation of its domain from

L_{0}

L_{1}

,while

η = \frac{[2 x - (L_{1} + L_{0})]}{L_{1} - L_{0}}

(51)

By substituting equation (49) into equations (29)–(35), and (46)–(48), the Lagrange energy functional

L_{Π}

associated with the CNT-reinforced laminated cylindrical shell is formulated as follows

L_{Π} = T_{c} + T_{p} - U_{c} - U_{p} - U_{s p r i n g} + W_{Δ T} + W_{Δ C} + W_{Δ P} + U_{g} + W_{E}

(52)

By applying the Ritz minimization procedure to equation (52) with respect to the unknown coefficients, the following result can be obtained

\frac{\partial L_{Π}}{\partial U_{m}} = \frac{\partial L_{Π}}{\partial V_{m}} = \frac{\partial L_{Π}}{\partial W_{m}} = 0

(53)

The discrete form of the governing equations can be written as follows

M \ddot{X} + C_{Δ p} \dot{X} + (K_{c} + K_{s p r i n g} + K_{Δ p} + K_{Δ T} + K_{Δ C} + K_{g} + K_{E}) X = 0

(54)

Solution of the governing equations

In the aeroelastic analysis, the governing equations are transformed into eigenvalue equations according to the frequency domain analysis method as

[Ω^{2} M + Ω C_{Δ p} \dot{X} + (K_{c} + K_{s p r i n g} + K_{Δ p} + K_{Δ T} + K_{Δ C} + K_{g} + K_{E})] X = 0

(55)

The prerequisite for the preceding equation to have a non-zero solution is

| Ω^{2} M + Ω C_{Δ p} + (K_{c} + K_{s p r i n g} + K_{Δ p} + K_{Δ T} + K_{Δ C} + K_{g} + K_{E}) | = 0

(56)

The complex eigenvalue $Ω = Ω_{R} + i ω$ is derived by solving equation (56). The eigenfrequency and the real part of the eigenvalue are named by $ω$ and $Ω_{R}$ , correspondingly. The determination of the system’s damping ratio can be achieved through the following expression

ξ = - \frac{Ω_{R}}{\sqrt{ω^{2} + Ω_{R}}}

(57)

It is widely recognized (Song and Li, 2012) that the point at which the damping ratio transitions from positive to negative signifies the critical flutter pressure $p_{\infty}^{c r}$ .

In the dynamic stability study, the cylindrical shells is considered to be placed under a periodically varying axial load $N_{x} (t)$ , and $N_{x} (t) = N_{c r} (u_{s} + u_{d} \cos ω_{d} t)$ , where $u_{s}$ , $u_{d}$ are named as the static and dynamic segments of load, correspondingly, and $ω_{d}$ is named as the loading frequency. The critical load for buckling of the shells is named as $N_{c r} = E h^{2} R \sqrt{3 - 3 v^{2}}$ (Qinkai and Fulei, 2013).

Substituting $N_{x} (t)$ into equation (54) and using the coordinate transformation $X = e^{(- C_{Δ p} / 2 M) t} δ$ to obtain the classical Mathieu-Hill equation as

M \ddot{δ} + C_{Δ p} \dot{δ} + (K_{c} + K_{s p r i n g} + K_{Δ p} + K_{Δ T} + K_{Δ C} + K_{E} + u_{s} N_{c r} K_{g} - C_{Δ p}^{2} / 4 M) δ + (u_{d} N_{c r} K_{g} \cos ω_{d} t) δ = 0

(58)

Equation (58) reveals the dynamically stable behavior of the cylindrical shells under periodic load, and its stable and unstable regions can be solved by Bolotin et al. (1965). For the symmetric cylindrical shells, instability region analysis reveals larger $2 T$ period domains than $T$ period counterparts, with the $2 T$ period solution formulated as

δ = a_{1} \sin \frac{k ω_{d} t}{2} + b_{1} \cos \frac{k ω_{d} t}{2}

(59)

Substituting equation (59) into equation (58) and equating respective sine/cosine coefficients of

a_{1}

and

b_{1}

yields their discrete linear algebraic system, requiring the coefficient matrix determinant to vanish. Solving the following equation determines the unstable region’s two boundary lines, whose intersection defines its origin

\det [(\begin{array}{c} K - \frac{1}{2} Q & 0 \\ 0 & K + \frac{1}{2} Q \end{array}) - ω_{d}^{2} (\begin{array}{c} \frac{1}{4} M & 0 \\ 0 & \frac{1}{4} M \end{array})] = 0

(60)

where

K = K_{c} + K_{s p r i n g} + K_{Δ p} + K_{Δ T} + K_{Δ C} + K_{E} + u_{s} N_{c r} K_{g} - C_{Δ p}^{2} / 4 M Q = u_{d} N_{c r} K_{g}

(61)

Parametric studies

This paper examines aeroelastic/dynamic instability in F-F multilayer cylindrical shells with CNT reinforcement under hygrothermal conditions. Where the CNT mass density is $ρ_{r} = 1400 k g / m^{3}$ , and Hill’s modulus of elasticity is $k_{r} = 271 G P a$ , $l_{r} = 88 G P a$ , $m_{r} = 17 G P a$ , $n_{r} = 1089 G P a$ , and $p_{r} = 442 G P a$ , respectively. For the reinforced matrix phase the material parameters are $ρ_{m} = 1150 k g / m^{3}$ , $β_{m} = 0.06$ , $α_{m} = (1 + 0.0005 Δ T) \times 1 0^{- 6} K$ and $v_{m} = 0.34$ . For the piezoelectric film the material parameters are $E_{11}^{p} = 30.34 G P a$ $E_{22}^{p} = 15.86 G P a$ , $G_{12}^{p} = 6.31 G P a$ , $ρ_{p} = 5516 k g / m^{3}$ , $ν_{12}^{p} = 0.31 Ξ_{33} = 1.43 \times 1 0^{- 8} F / m$ , $d_{31} = - 170 \times 1 0^{- 12} m / V$ , $d_{32} = - 100 \times 1 0^{- 12} m / V$ .

Equation validation

In this section, the derived formulations and computational algorithms were validated through benchmark comparisons against previous literature cases. This verification confirms both analytical correctness and computational efficacy, to provide a rigorous foundation for subsequent parametric investigation of dynamic characteristics. In terms of the first benchmark, the material parameters of the composite cylindrical shell are

E_{1} = 25 E_{2}

E_{2} = 2 G P a

G_{12} = 0.5 E_{2}

v_{12} = 0.25

and

ρ = 1500 k g / m^{3}

, while the dimensionless frequency parameter of the shells can be expressed as

Ω = ω L^{2} / h \sqrt{ρ / E_{2}}

. The comparison dates of the dimensionless natural frequency are detailed in Table 1, where three case of boundary conditions are studied, F-C indicates the shells is free at

x = 0

and clamped at

x = L

, F-SS is left end of free and right end of simply supported, and C-SS is left clamped and right simply supported. The spring stiffness values corresponding to the different conditions of boundary are, F:

k_{u} = k_{v} = k_{w} = K_{w} = 0

; SS:

k_{u} = k_{v} = k_{w} = 1 0^{4} E_{2}

K_{w} = 0

; SD:

k_{u} = 0

k_{v} = k_{w} = 1 0^{4} E_{2}

K_{w} = 0

and C:

k_{u} = k_{v} = k_{w} = K_{w} = 1 0^{4} E_{2}

. SD indicates boundary condition of soft simply supported. For the second instance, the properties of material of the composite cylindrical shells are

E_{1} = 20 E_{2}

E_{2} = 10 G P a

G_{12} = 0.65 E_{2}

v_{12} = 0.25

, and

ρ = 1500 k g / m^{3}

, while the dimensionless frequency parameter of the shells with cylindrical geometry is

Ω = ω R \sqrt{ρ / E_{2}}

Table 1.

Dimensionless natural frequency of the shells with cylindrical geometry. ( $L = 5 m$ , $R = 1 m$ ; $h = 0.05 m$ ; $[0^{\circ} / 9 0^{\circ} / 0^{\circ}]$ ).

$m (n = 1)$	F-C			F-SS			C-SS
$m (n = 1)$	(Li et al., 2019a)	Present	Error	(Li et al., 2019a)	Present	Error	(Li et al., 2019a)	Present	Error
$1$	$74.347$	$74.583$	$0.3 %$	$73.605$	$73.436$	$- 0.2 %$	$157.652$	$159.355$	$1.1 %$
$2$	$224.296$	$224.296$	$2.5 %$	$221.642$	$226.464$	$2.2 %$	$314.168$	$318.263$	$1.3 %$
$3$	$389.134$	$395.542$	$1.6 %$	$384.699$	$385.509$	$1.3 %$	$472.285$	$476.689$	$0.9 %$
$4$	$542.657$	$552.548$	$1.8 %$	$536.699$	$543.691$	$1.3 %$	$628.262$	$635.453$	$1.1 %$

The numerical results in Table 1 and Table 2 are in good agreement with the literature results in both the present derivation methodology and the computational reliability of the implemented algorithms.

Table 2.

Dimensionless natural frequency of the cylindrical shells ( $L = 4 m$ ; $R = 1 m$ ; $h = 0.01 m$ ; $[1 5^{\circ} / 1 5^{\circ} / 1 5^{\circ} / 1 5^{\circ}]$ ).

	$m = 4, n = 3$
	F-F	C-C	SD−SD
Qu et al. (2013)	$0.69183$	$0.86438$	$0.85660$
Present	$0.69115$	$0.86587$	$0.85763$
Error	$- 0.1 %$	$0.1 %$	$0.1 %$

The comparison results of the natural frequency of double-ended clamped cylindrical shells in a hygrothermal environment is presented in Table 3. The environmental parameters are set as

T = 350 K

C = 0.5 %

. The geometric parameters of the cylindrical shell are defined as

L / R = 5

h / R = 0.01

. The material properties include

E_{1} = 130 G P a

E_{2} = 8 G P a

G_{12} = 5.5 G P a

ν_{21} = 0.3

and

ρ = 1600 k g / m^{3}

. Additionally, the thermal expansion coefficients are

α_{1} = - 0.3 \times 1 0^{- 6} /^{\circ} C

and

α_{2} = 28.1 \times 1 0^{- 6} /^{\circ} C

, while the moisture expansion coefficients are

β_{1} = 0

and

β_{2} = 0.44 \times 1 0^{- 6}

. Comparisons show this study’s results align with literature, validating the derived formulations and temperature- and humidity-related programming.

Table 3.

Natural frequency (Hz) of cylindrical shells in a hygrothermal environment ( $T = 350 K$ ; $Δ C = 0.5 %$ ; $L / R = 5$ ; $h / R = 0.01$ ; $[0^{\circ} / 9 0^{\circ} / 0^{\circ}]$ ).

$m (n = 1)$	(Li et al., 2019b)		Present
$m (n = 1)$	Galerkin method	FEM	Present
$1$	$290.74$	$287.16$	$285.88$
$2$	$127.48$	$124.99$	$126.79$
$3$	$62.58$	$60.37$	$62.14$

Subsequently, the verification results for the flutter boundary are summarized in Table 4. The geometric parameters of the cylindrical shell used in this case are as follows:

L = 0.5 m

R = L / 4

h = R / 1000

. The material properties include

ρ = 1643 k g / m^{3}

μ_{12} = 0.3

E_{1} = 150 G P a

E_{2} = 9 G P a

G_{12} = 7.1 G P a

. For the dynamic instability verification, the results are presented in Table 5, where the geometric configuration is defined by

L / R = 2

and

R / h = 200

, and the material properties are given as

E_{1} / E_{2} = 40

G_{12} / E_{2} = 0.5

and

ν_{12} = 0.25

. As shown in tables, the flutter boundaries and onset points of dynamic instability calculated in this study are in close agreement with those reports in the literature.

Table 4.

Flutter boundaries for circular cylindrical shells with simply supported boundary conditions ( $L = 0.5 m$ ; $R = L / 4$ ; $h = R / 1000$ ).

	S-S
Lamination scheme	${[4 5^{\circ} / - 6 0^{\circ} / 6 0^{\circ}]}_{s}$	${[6 0^{\circ} / 4 5^{\circ} / 6 0^{\circ}]}_{s}$
(Chai et al., 2018)	$3151$	$3067$
Present	$3175$	$3038$
Error	$0.8 %$	$- 0.9 %$

Table 5.

Instability origins ( $ω^{*} = ω_{d} 2 π R \sqrt{ρ / A_{11}}$ ) in simply supported circular cylindrical shells ( $L / R = 2$ ; $R / h = 200$ ; $u_{s} = 0.1$ ).

Lamination scheme	Point of origin $ω^{*}$
Lamination scheme	(Ng et al., 1998)	(Yadav et al., 2022)	Present
$[0^{\circ} / 9 0^{\circ}]$	$0.5729$	$0.5603$	$0.5637$
$[9 0^{\circ} / 0^{\circ}]$	$0.5654$	$0.5681$	$0.5659$
$[9 0^{\circ} / 0^{\circ} / 9 0^{\circ} / 0^{\circ}]$	$0.6714$	$0.6702$	$0.6756$

Aeroelasticity research

The aeroelasticity of laminated cylindrical shells was investigated with the initial volume fraction of CNT reinforced fibers is $f_{r} = 0.1$ , the geometrical size of the cylindrical shells are $L = 1 m$ , $R = 1 m$ , and $h = 0.005 m$ ; the initial airflow related parameters are $γ = 1.4$ , $a_{\infty} = 213.36 m / s$ , $M_{a} = 2$ , and $Λ = 0^{\circ}$ , the initial temperature and humidity are $300 K$ and $0.1 %$ , the layup of the cylindrical shells are $3$ layers, and the layup angle is $[Θ^{\circ} / Θ^{\circ} / Θ^{\circ}]$ .

This study first investigated the convergence behavior for flutter boundaries of cylindrical shells under supersonic flow. Figure 2 shows the evolution of critical flutter points versus axial modal truncation order. It can be seen that when the axial truncation order reaches $12$ and $13$ , the change of the zero-damping point is very small, which shows a good stable convergence result. Subsequently, it is reasonable to select the axial truncation order of $13$ .

Figure 2.

Flutter boundary under varying axial truncation orders.

To quantitatively elucidate the impact of carbon nanotube volume fraction, the critical flutter pressure corresponding to various volume fractions is presented in Figure 3. In this figure, we take $f_{r} = 0.10, 0.12, 0.14$ and $0.16$ . An escalation in the volume fraction results in a corresponding rise in the critical aerodynamic pressure. The reason is increasing the CNT volume fraction might strengthen flexural rigidity through more effective reinforcement mechanisms. The resultant stiffening of the laminated cylindrical shell consequently expands the stability range.

Figure 3.

Dependence of flutter criticality on carbon nanotube volume fraction ( $Θ = 4 5^{\circ}$ ).

The flutter pressure corresponding to the angle of carbon nanotube fibers is presented in Figure 4. It is seen that the critical point of flutter exhibits a gradual intensify with the layup angle decreasing. This phenomenon is predominantly governed by supersonic axial flow dynamics. Reducing the carbon nanotube fiber orientation angle enhances axial structural rigidity in laminated cylindrical shells. Consequently, CNT cylindrical shells with smaller ply angles demonstrate superior aeroelastic stability characteristics.

Figure 4.

Dependence of flutter criticality on $Θ$ .

Figure 5 presents the parametric dependencies of critical flutter pressure in composite cylindrical shells, through systematic examination of airflow deflection angle variations and Mach number effect. For example, a $Θ = 4 5^{\circ}$ $f_{r} = 0.1$ cylindrical shell is taken into consideration. Four sets of Mach numbers of $2.0 M_{a}$ , $2.5 M_{a} 3.0 M_{a}$ and $3.5 M_{a}$ are involved. It is observed that for all Mach numbers, the critical aerodynamic pressures increase as the airflow deflection angle increases. According to Figure 1 and the influence of CNT ply orientation on critical dynamic pressure, it shows that at $4 5^{\circ}$ CNT ply orientation, increased airflow shell axis angle effectively decreases the angular deviation between airflow and reinforcement fibers. This improved fiber load alignment can improve CNT reinforcement efficacy and structural stiffness. Conversely, higher Mach numbers reduce the flutter pressure at any airflow deflection angle.

Figure 5.

Dependence of flutter criticality on $Λ$ ( $Θ = 4 5^{\circ}$ , $f_{r}$ ).

The impact of hydrothermal environment on the critical aerodynamic pressures of the hierarchical shells with cylindrical geometry are illustrated in Figure 6. The findings demonstrate reduced flutter boundaries for hierarchical cylindrical shells under hygrothermal exposure. This degradation originates from temperature- and moisture-induced reduction in elastic modulus, which compromises structural stiffness and consequently destabilizes the system.

Figure 6.

Dependence of flutter criticality on the hydrothermal environment ( $Θ = 4 5^{\circ}$ ).

Figure 7 quantifies the relationship between critical aerodynamic pressure and piezoelectric thickness ratio $(h_{p} / h)$ in cylindrical hierarchical shells with surface-bonded piezoelectric films. The critical pressure exhibits direct proportionality to $h_{p} / h$ , with progressive increases in thickness ratio yielding exponentially augmented pressure resistance. The main reason is that the increase in macro-fiber composite (MFC) thickness significantly improves the bending stiffness of the hierarchical cylindrical shell

Figure 7.

Effect of thickness ratio on critical flutter dynamic pressure.

Figure 8 demonstrates that MFC actuator angular orientation significantly affects cylindrical shells’ critical aerodynamic pressure. At $0.02$ thickness ratio and $300 V$ voltage, the flutter boundary increases with piezoelectric fiber angle variation, peaking at $α = 9 0^{\circ}$ .

Figure 8.

Effect of MFC angle on critical aerodynamic pressure.

Figure 9 depicts the influence of electric field voltage on flutter pressure for varying MFC fiber orientations in cylindrical hierarchical shells at thickness ratio of $0.02$ . The data demonstrate that an initial increase followed by progressive decay within the range of $- 30 V$ to $360 V$ . Especially, equivalent critical pressures occur at $0 V$ bias regardless of fiber orientation, which is a baseline condition matching the literature (Zhou et al., 2019). This voltage-neutral response confirms angle independence at zero potential.

Figure 9.

Critical aerodynamic pressure as variation of voltage and MFC angle.

Dynamic instability

As discussed in section 2.6 for cylindrical shells subjected to periodic parametric loading, the unstable regions of the hierarchical cylindrical shells coupled with MFC actuators exhibit a complex dependence on multiple parameters. These include aerodynamic loads, Mach number, airflow deflection angle, CNT layup angle, volume fraction, temperature, humidity, static load and applied electric field voltage. The dimensionless critical excitation frequency is named as $ω^{*} = ω_{d} L^{2} / h \sqrt{ρ / E_{2}}$ .

Figure 10(a) illustrates aerodynamic pressure effects on instability regions in cylindrical hierarchical shells under F-F boundary conditions. The V-shaped area formed by three curves delineates unstable regions across varying conditions, with exterior zones indicating stable regions. It can be seen that variations in aerodynamic pressure substantially affect both the extent and positioning of unstable regions. Specifically, an increase in aerodynamic pressure leads to a notable shift of the regions of instability in the direction of the high-frequency band. Elevated aerodynamic pressure increases the natural frequency of laminated composite structures within their operational range, which consequently raise the critical transverse excitation frequency. Figure 10(b) illustrates the influences of three different Mach numbers on the regions of instability of the hierarchical shells. The variation shown in Figure 10 indicates a negative correlation between the shift in the instability domain and the rise in Mach number.

Figure 10.

Effects of aerodynamic pressure and Mach number on the regions of instability. (a) $Θ = 4 5^{\circ}, M_{a} = 2.0$ (b) $Θ = 4 5^{\circ}, p_{\infty} = 100 k P a$ .

Figure 11(a) demonstrates that the airflow deflection angle exerts negligible influence on instability boundaries. However, the ply angle of the hierarchical shells with cylindrical geometry has a substantial influence on the regions of instability, as shown in Figure 11(b). The results demonstrate that decreasing the ply angle simultaneously might change instability regions toward higher frequencies while reducing their bandwidth. This phenomenon is primarily attributed to the layup configuration, which significantly influences the structural stiffness. In this study, the reduction in ply orientation angle leads to an increase in the structure’s natural frequency, consequently elevating the corresponding critical excitation frequency.

Figure 11.

Effects of airflow angle and ply-angle on the regions of instability. (a) $Θ = 4 5^{\circ}, p_{\infty} = 100 k P a, M_{a} = 2.0$ . (b) $Λ = 0^{\circ}, p_{\infty} = 100 k P a, M_{a} = 2.0$ .

Figure 12(a) examines how the percentage of CNTs influences the unstable regions of the hierarchical shells. As illustrated in the figure, the volume fraction has a substantial impact on the extent of the regions of instability. As the volume fraction increases, the regions of instability shift to a high-frequency band. This phenomenon can be attributed to the fact that the increase in volume fraction improves the bending stiffness of the system. However, the width of this region of instability remains approximately constant. Figure 12(b) demonstrates static load parameter effects on instability regions in cylindrical hierarchical shells. Increased static loads induce backward shifts of instability regions, with concomitant reductions in V-shaped enclosed unstable areas.

Figure 12.

Variation of the regions of instability with volume fraction and static load. (a) $Θ = 4 5^{\circ}, p_{\infty} = 100 k P a, u_{s} = 0.1$ (b) $Θ = 4 5^{\circ}, p_{\infty} = 100 k P a, f_{r} = 0.1$ .

Figure 13 displays instability regions under four temperature/humidity conditions, showing that maximum instability corresponds to elevated temperature and humidity levels, occurring predominantly at lower frequency. This is due to the reduction in the system’s flexural rigidity resulting from the increase in temperature and humidity.

Figure 13.

Dependence of unstable regions on temperature and humidity. (a) $Θ = 4 5^{\circ}, p_{\infty} = 100 k P a, C = 0.1 %$ . (b) $Θ = 4 5^{\circ}, p_{\infty} = 100 k P a, T = 300 K$ .

To further illustrate the effect of MFC on the structural regions of instability, Figure 14 demonstrates that increasing the fiber orientation angle ( $0^{\circ} - 9 0^{\circ}$ ) in MFC-integrated hierarchical shells simultaneously influences instability regions toward higher frequencies while reducing their spatial extent, with both trends quantified parametrically across the angular range.

Figure 14.

Dependence of unstable regions on MFC angle. (a) and (b) respectively represent the three-dimensional and two-dimensional plane diagrams of the dynamic instability region varying with the MFC angle. ( $Θ = 4 5^{\circ}, p_{\infty} = 100 k P a, T = 300 K, C = 0.1 %, V_{0} = 300 V$ ).

Figure 15 illustrates the implications of voltage on the regions of instability of structural when parameters $α = 4 5^{\circ}$ and $h_{p} / h = 0.02$ . As illustrated in Figure 15(a), for negative voltages, instability domains migrate toward higher frequencies while progressively contracting as voltage magnitude increases. However, positive voltages drive instability toward lower frequencies while expanding its spatial footprint. Transverse polarization in piezoelectric fibers causes the polarity-dependent behavior. Negative voltages induce transverse compression that reduces bending stiffness, while positive voltages generate transverse tension that strengthens structural rigidity. This phenomenon is consistent with findings in the literature (Eftekhar et al., 2018).

Figure 15.

Dependence of unstable regions on voltage. (a) and (b) respectively represent the three-dimensional and two-dimensional plane diagrams of the dynamic instability region varying with voltage. ( $Θ = 4 5^{\circ}, p_{\infty} = 100 k P a, T = 300 K, C = 0.1 %, α = 4 5^{\circ}$ ).

Conclusion

This study aims to investigate the aeroelasticity and dynamic instability of multi-layered cylindrical shells reinforced with CNT composites under F-F boundary conditions. We quantify the effects of airflow parameters, CNT volume fraction, fiber ply orientation, temperature, humidity, and voltage on the flutter boundary and dynamic instability boundary. Key findings can be summarized as below.

1. Increasing CNT fiber volume fraction, MFC fiber angle, and thickness ratio while decreasing CNT fiber ply angle, might elevate the flutter and dynamic instability boundaries of double-ended free cylindrical shells, thereby improving the stability of elastic structures. These findings provide valuable insights for the structural design of high-speed flight projectiles from a vibration perspective.

2. Elevated temperature, humidity, and Mach number compromise the flutter and dynamic instability thresholds of double-ended free cylindrical shells. While increased aerodynamic pressure can enhance dynamic instability resistance, this benefit is counterbalanced by its tendency to promote flutter onset. Consequently, missile structural design requires careful multi-parameter optimization to mitigate these competing aeroelastic effects.

3. Voltage variations significantly influence on both flutter and dynamic instability boundaries in double-ended free cylindrical shells. This voltage-dependent control mechanism directly provides future active vibration controller designs for such aerospace structures.

Footnotes

ORCID iD

Yanhong Wu

Funding

The authors disclosed receipt of the following financial support for the research, authorship, and/or publication of this article: The project was supported by the key R&D plan of the province of Heilongjiang. Project No: 2022ZX01A14.

Declaration of conflicting interests

The authors declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.

References

Adams

Miller

(1977) Hygrothermal microstresses in a unidirectional composite exhibiting inelastic material behavior. Journal of Composite Materials 11: 285–299.

Arani

Kiani

Afshari

(2019) Aeroelastic analysis of laminated FG-CNTRC cylindrical panels under yawed supersonic flow. International Journal of Applied Mechanics 11(06): 1950052.

Asadi

Wang

(2017) Dynamic stability analysis of a pressurized FG-CNTRC cylindrical shell interacting with supersonic airflow. Composites Part B: Engineering 118: 15–25.

Avramov

Chernobryvko

Uspensky

, et al. (2019) Self-sustained vibrations of functionally graded carbon nanotubes-reinforced composite cylindrical shells in supersonic flow. Nonlinear Dynamics 98(3): 1853–1876.

Bisheh

Rabczuk

(2020) Free vibration analysis of smart laminated carbon nanotube-reinforced composite cylindrical shells with various boundary conditions in hygrothermal environments. Thin-Walled Structures 149: 106500.

Bolotin

Weingarten

Greszczuk

, et al. (1965) Dynamic stability of elastic systems. Journal of Applied Mechanics 32: 718.

Boukhoulda

Adda-Bedia

Madani

(2006) The effect of fiber orientation angle in composite materials on moisture absorption and material degradation after hygrothermal ageing. Composite Structures 74(4): 406–418.

Bowles

Tompkins

(1989) Prediction of coefficients of thermal expansion for unidirectional composites. Journal of Composite Materials 23(4): 370–388.

Chai

Y-Y

Song

Z-G

F-M

(2017) Active aerothermoelastic flutter suppression of composite laminated panels with time-dependent boundaries. Composite Structures 179: 61–76.

10.

Chai

Song

(2018) Investigations on the aerothermoelastic properties of composite laminated cylindrical shells with elastic boundaries in supersonic airflow based on the rayleigh–ritz method. Aerospace Science and Technology 82-83: 534–544.

11.

Chakraborty

Dey

Kumar

(2021) Instability characteristics of damped CNT reinforced laminated shell panels subjected to in-plane excitations and thermal loading. Structures 34: 2936–2949.

12.

Dat

Quan

Mahesh

, et al. (2020) Analytical solutions for nonlinear magneto-electro-elastic vibration of smart sandwich plate with carbon nanotube reinforced nanocomposite core in hygrothermal environment. International Journal of Mechanical Sciences 186: 105906.

13.

Eftekhar

Zeynali

Nasihatgozar

(2018) Electro-magneto temperature-dependent vibration analysis of functionally graded-carbon nanotube-reinforced piezoelectric mindlin cylindrical shells resting on a temperature-dependent, orthotropic elastic medium. Mechanics of Advanced Materials and Structures 25(1): 1–14.

14.

Hasheminejad

Aghayi Motaaleghi

(2015) Aeroelastic analysis and active flutter suppression of an electro-rheological sandwich cylindrical panel under yawed supersonic flow. Aerospace Science and Technology 42: 118–127.

15.

Heidari-Soureshjani

Rajabi

Talebitooti

, et al. (2024a) Hygrothermal vibro-buckling of FG ceramic-steel porous consolidated conical-conical shells. Thin-Walled Structures 201: 112002.

16.

Heidari-Soureshjani

Roohollah

Talebitooti

(2024b) On the frequency characteristics of rotating combined conical-conical shells made of FG-CNTRC composite materials under thermal environments. Mechanics Based Design of Structures and Machines 52(3): 1437-1461.

17.

Heydarpour

Malekzadeh

(2018) Dynamic stability of rotating FG-CNTRC cylindrical shells under combined static and periodic axial loads. International Journal of Structural Stability and Dynamics 18(12): 1850151.

18.

Hill

(1965) A self-consistent mechanics of composite materials. Journal of the Mechanics and Physics of Solids 13(4): 213–222.

19.

(2018) Parametric instability of a functionally graded cylindrical thin shell subjected to both axial disturbance and thermal environment. Thin-Walled Structures 123: 25–35.

20.

Pang

Miao

, et al. (2019a) A semi analytical method for free vibration analysis of composite laminated cylindrical and spherical shells with complex boundary conditions. Thin-Walled Structures 136: 200–220.

21.

Xie

(2019b) Vibration characteristics of a rotating composite laminated cylindrical shell in subsonic air flow and hygrothermal environment. International Journal of Mechanical Sciences 150: 356–368.

22.

Mohammadi

Ashenai Ghasemi

Shahgholi

(2021) Stability analysis of an axially moving nanocomposite circular cylindrical shell with time-dependent velocity in thermal environments. Mechanics Based Design of Structures and Machines 49(5): 659–688.

23.

Lam

Reddy

(1998) Dynamic stability of cross-ply laminated composite cylindrical shells. International Journal of Mechanical Sciences 40(8): 805–823.

24.

Pakravan

Soureshjani

Talebitooti

, et al. (2022) Haar wavelet technique applied on the functionally graded carbon nanotube reinforced conical shells to study free vibration and buckling behaviors in thermal environments. Journal of Vibration and Control 28(15-16): 1863–1878.

25.

Qinkai

Fulei

(2013) Parametric instability of a rotating truncated conical shell subjected to periodic axial loads. Mechanics Research Communications 53: 63–74.

26.

Long

, et al. (2013) A unified formulation for vibration analysis of composite laminated shells of revolution including shear deformation and rotary inertia. Composite Structures 98: 169–191.

27.

Song

Z-G

F-M

(2012) Active aeroelastic flutter analysis and vibration control of supersonic composite laminated plate. Composite Structures 94(2): 702–713.

28.

Wang

Song

Zang

, et al. (2023) Experimental and theoretical investigation on vibration of laminated composite conical-cylindrical-combining shells with elastic foundation in hygrothermal environment. Composite Structures 323: 117470.

29.

Yadav

Amabili

Panda

, et al. (2022) A semi-analytical approach for instability analysis of composite cylindrical shells subjected to harmonic axial loading. Composite Structures 296: 115882.

30.

Zhang

Song

Liew

(2017) Modeling aerothermoelastic properties and active flutter control of nanocomposite cylindrical shells in supersonic airflow under thermal environments. Computer Methods in Applied Mechanics and Engineering 325: 416–433.

31.

Zhou

Yang

(2019) Aeroelastic stability analysis of curved composite panels with embedded macro fiber composite actuators. Composite Structures 208: 725–734.