曲线纤维变刚度壁板的亚音速气弹稳定性研究

段静波; 徐步青

doi:10.3969/j.issn.0258-2724.20200277

曲线纤维变刚度壁板的亚音速气弹稳定性研究

doi: 10.3969/j.issn.0258-2724.20200277

段静波,
徐步青^,

石家庄铁道大学工程力学系，河北石家庄 050043

基金项目: 国家自然科学基金(11702325)；河北省重点研发计划（21350401D）

详细信息

作者简介:
段静波（1982—），男，副教授，研究方向为气动弹性力学，E-mail：duanjingbo@stdu.edu.cn

通讯作者:
徐步青（1970—），男，副教授，研究方向为气动弹性力学，E-mail：623220972@126.com

中图分类号: V221.3
计量
- 文章访问数: 393
- HTML全文浏览量: 205
- PDF下载量: 17
- 被引次数: 0
出版历程
- 收稿日期: 2020-05-18
- 修回日期: 2021-02-06
- 网络出版日期: 2022-08-18
- 刊出日期: 2021-03-03

Aeroelastic Instability of Variable-Stiffness Panels with Curvilinear Fibers in Subsonic Flow

DUAN Jingbo,
XU Buqing^,

Department of Engineering Mechanics, Shijiazhuang Tiedao University, Shijiazhuan 050043, China

摘要

摘要:
针对曲线纤维变刚度复合材料层合板在高速列车壁板结构轻量化设计中的广阔应用前景，研究了弹性、黏弹性变刚度复合壁板在亚音速流场中的气动弹性稳定性问题. 首先，基于Mindlin厚板理论和势流流动理论分别描述壁板结构变形和亚音速气动力，根据虚功原理和有限元法建立了曲线纤维变刚度复合材料弹性/黏弹性壁板的气动弹性稳定性分析模型，进而采用复模态理论求解复合材料变刚度壁板发散临界速度；在验证方法正确性和收敛性的基础上，研究了壁板关键参数等对复合变刚度壁板发散失稳特性的影响规律. 研究结果表明：与直线纤维壁板相比，通过调整曲线纤维路径可以实现壁板的发散临界速度50%左右的提升，有效增强壁板气动弹性稳定性.
- 气动弹性稳定性 /
- 曲线纤维 /
- 变刚度 /
- 复合材料壁板 /
- 亚音速气流
Abstract:
In view of the extensive application of curved fiber composite laminates in the lightweight design of high-speed train structures, the aeroelastic stability of elastic and viscoelastic variable-stiffness composite panels in a subsonic flow field was studied. First, classical thick theory along with a Mindlin plate was adopted for structural modeling and potential flow theory for aerodynamic modeling. An aeroelastic model of composite variable-stiffness panels with curvilinear fibers was then established adopting the principle of virtual work and the finite element method, which was solved using complex mode theory in the frequency domain. The divergence characteristics for key parameters were investigated following verification of the validity and convergence of the presented method. Numerical results show that, relative to the straight-fiber panel, the critical divergence speed can be increased by approximately 50% by varying the path orientations of the curvilinear fibers.
- aeroelastic stability /
- curvilinear fibers /
- variable stiffness /
- composite panels /
- subsonic flow

HTML全文

为追求列车设计的高速化与轻量化，新型复合材料越来越广泛地应用于高速轨道交通领域，复合材料曲线纤维层合板作为一种近年新出现的复合材料，与传统的直线纤维铺层相比，由于纤维方向随位置的变化而变化，曲线纤维铺层呈现出变刚度的特点. Gűrdal等^[1]研究表明，刚度在空间上非均匀分布的变刚度壁板能大幅改善结构承载能力，适应更为复杂的力学环境. 将曲线纤维变刚度复合材料层合板应用于高速列车蒙皮、裙板等壁板结构，有利于新一代复合材料高铁列车高速化与轻量化设计目标的顺利实现.

然而，伴随着高速列车行进速度的不断提高，车体蒙皮、裙板等壁板类结构所面临的力学环境越来越恶劣，设计不合理可能会导致壁板发生剧烈的气动弹性振动，严重影响结构的使用寿命甚至列车的运行安全. 因此，亚音速曲线纤维变刚度壁板气动弹性稳定性问题将是车体新型复合材料应用中的一个关键问题.

曲线纤维复合材料最早由Hyper等^[2]提出，Gűrdal等^[1,3]则明确提出了变刚度复合材料的概念. 随后，学者们对于曲线纤维变刚度壁板力学问题进行广泛研究：Hamed等^[4]对曲线纤维变刚度壁板进行大变形和应力分析；马洪涛^[5]开展了曲线纤维变刚度壁板的力学性能研究；聂国隽等^[6]建立了曲线纤维增强复合材料层合板自由振动问题的基本方程，得到了层合板的自振频率及相应的振型；马成^[7]通过对纤维曲线铺放层合板阻尼比的研究，探究纤维曲线铺放层合板的减振效果；Groh等^[8]研究了横向剪切效应对变刚度层合板屈曲性能的影响；Hao等^[9]研究了复合材料变刚度板的屈曲及优化问题，提出了复合材料变刚度结构的精确建模、分析与结构优化一体化设计框架；孙士平等^[10]研究了复合载荷作用下变刚度复合材料回转壳的屈曲优化问题.

迄今为止，有关曲线纤维变刚度壁板气动弹性问题的研究较少，但已经引起了国内外学者的关注：Vahid等^[11]研究了高速流场中的曲线纤维变刚度层合壁板颤振问题；欧阳小穗和刘毅^[12]研究了高速流场中的曲线纤维变刚度层合壁板非线性颤振响应；Touraj等^[13]将曲线纤维变刚度壁板应用于飞机机翼承力结构，研究了机翼-引擎系统的气动弹性特性.

基于此，本文以新一代新型复合材料高速列车壁板结构为背景，基于目前列车设计时速范围，研究曲线纤维变刚度壁板的亚音速气动弹性稳定性问题，考虑壁板可能出现的黏弹性效应，分析复合材料关键参数对复合变刚度壁板气动弹性稳定性的影响规律.

1. 曲线纤维路径表征

图1为曲线纤维复合材料层合板示意，以层合板中第k层单层变刚度铺层为研究对象(图1（b）)，图中：${{\textit{z}}_k}$和${{\textit{z}}_{k + 1}}$分别为该铺层上、下表面在坐标系xOz中的坐标值；$ {h_k} $为该铺层厚度；$ {V_\infty } $为来流速度；L、H、h分别为矩形变刚度复合材料层合板的长度、宽度、厚度；在局部坐标系$ x'Oy' $下，T₀为纤维在板中心处与$ x' $轴方向的夹角；T₁为纤维在边界$\left( x' = \pm {L}/{2} \right)$处与$ x' $轴方向的夹角；ϑ为$ x' $轴与x轴的夹角，表示该铺层的方向角. 该铺层曲线纤维方向角用 <T₀｜T₁>表示. 假设纤维方向从板的中心位置开始呈线性变化，则该铺层在局部坐标系下任意位置处纤维方向角可表示为

图 1 曲线纤维层合板

Figure 1. Curvilinear fibers composite panels

下载: 全尺寸图片幻灯片

$$ T(x') = \frac{{2({T_1} - {T_0})}}{L}\left| {x'} \right| + {T_0} . $$

(1)

2. 复合变刚度壁板颤振方程

2.1 壁板变形几何关系

结构方面，考虑厚度方向上的剪切变形，忽略厚度方向上法向变形，采用Mindlin厚板理论描述曲线纤维变刚度层合壁板，其应变-位移关系为

$$ \begin{gathered} {\boldsymbol{\varepsilon }} = \left[ {\begin{array}{*{20}{c}} {{{\boldsymbol{\varepsilon }}_{{0}}}} \\ {\boldsymbol{\gamma }} \end{array}} \right] + \left[ {\begin{array}{*{20}{c}} {{\textit{z}}{\boldsymbol{\kappa }}} \\ {\mathbf{0}} \end{array}} \right] , \\ \end{gathered} $$

(2)

式中：$ {\boldsymbol{\varepsilon}} $为线应变；$ {{\boldsymbol{\varepsilon }}_{\text{0}}} $为壁板中面面内位移产生的应变向量，如式（3）；$ {\boldsymbol{\gamma }} $为横向剪切应变向量，如式（3）；$ {\boldsymbol{\kappa }} $为弯曲时壁板的曲率向量，如式（3）；z为壁板厚度方向的坐标.

$$ \left\{ \begin{array}{l} {{\boldsymbol{\varepsilon }}_0} = {\left[ {\begin{array}{*{20}{l}} {\dfrac{{\partial u}}{{\partial x}}}&{\dfrac{{\partial \upsilon }}{{\partial y}}}&{\dfrac{{\partial u}}{{\partial y}} + \dfrac{{\partial \upsilon }}{{\partial x}}} \end{array}} \right]^{\rm{T}}},\\ {\boldsymbol{\kappa }} = {\left[ {\begin{array}{*{20}{c}} {\dfrac{{\partial {\theta _y}}}{{\partial x}}}&{\dfrac{{\partial {\theta _x}}}{{\partial y}}}&{\dfrac{{\partial {\theta _x}}}{{\partial x}} + \dfrac{{\partial {\theta _y}}}{{\partial y}}} \end{array}} \right]^{\rm{T}}},\\ {\boldsymbol{\gamma }} = {\left[ {\begin{array}{*{20}{c}} {\dfrac{{\partial w}}{{\partial x}} - {\theta _x}}&{\dfrac{{\partial w}}{{\partial y}} - {\theta _y}} \end{array}} \right]^{\rm{T}}}, \end{array} \right. $$

(3)

式中：$ u $、$ \upsilon $分别为中面上的点沿x轴和y轴方向的位移；w为z轴方向的位移；$ {\theta _x} $、$ {\theta _y} $分别为中面绕x轴和y轴方向的转角.

2.2 变刚度壁板弹性/黏弹性本构关系

对于正交各向异性铺层材料，采用kelvin-Voigt黏弹性本构模型，则其应力-应变关系^[14]可写为

$$ \left\{ \begin{gathered} \mathit{\bf{\sigma }} = {{\boldsymbol{Q}}_{\rm{m}}}\left({\boldsymbol{\varepsilon }} + \eta {\boldsymbol{\dot \varepsilon }}\right) , \\ {\mathbf{\tau }} = {{\boldsymbol{Q}}_{\rm{s}}}\left({\boldsymbol{\gamma }} + \eta {\boldsymbol{\dot \gamma }}\right) , \\ \end{gathered} \right. $$

(4)

式中：$ \eta $为材料黏弹性阻尼系数，当$ \eta = 0 $时，铺层材料就退化为弹性材料；$\bf\sigma 、\tau 、\gamma 、\dot{\varepsilon }$分别为正应力、剪应力、切应变、线应变率；$ {{\boldsymbol{Q}}_{\text{m}}}、{{\boldsymbol{Q}}_{\text{s}}} $分别为面内刚度和剪切刚度，其元素与弹性本构关系一致，在此不再赘述，如（5）.

$$\left\{ { \begin{array}{l} {{\boldsymbol{Q}}}_{{\text{m}}}=\left[\begin{array}{ccc}{Q}_{11}& {Q}_{12}& 0\\ {Q}_{12}& {Q}_{22}& 0\\ 0& 0& {Q}_{66}\end{array}\right]\text{，}\\ {{\boldsymbol{Q}}}_{\text{s}}=\left[\begin{array}{cc}{Q}_{44}& 0\\ 0& {Q}_{55}\end{array}\right]. \end{array} } \right.$$

(5)

由于曲线纤维角度影响，单层板内本构关系随壁板平面是变化的，当多个单层板形成层合板后，则曲线纤维角度和各层铺层角共同影响壁板本构关系，第k层本构关系为

$$ \left\{ \begin{array}{l} {\sigma }_{k}={\boldsymbol{\bar Q}}_{{\text{m}}k}({\varepsilon }_{k} + \eta {\dot{\varepsilon }}_{k})\text{，} \\ {\tau }_{k}={\boldsymbol{\bar Q}}_{{\text{s}}k}({\gamma }_{k} + \eta {\dot{\gamma }}_{k}) , \end{array} \right. $$

(6)

式中：$ {\sigma }_{k}、{\tau }_{k} $分别为该铺层的正应力和切应力；$ {\varepsilon }_{k}、{\gamma }_{k} $分别为该铺层的线应变和切应变；${\boldsymbol{ Q}}_{{\text{m}}k}、{\boldsymbol{ Q}}_{{\text{s}}k}$分别为偏轴面内刚度矩阵和剪切刚度矩阵，${{\boldsymbol{\bar Q}}_{{\text{m}}k}} = {{\boldsymbol{T}}_{{\text{m}}k}}{{\boldsymbol{Q}}_{\text{m}}}{\boldsymbol{T}}_{{\text{m}}k}^{\text{T}},{{\boldsymbol{\bar Q}}_{{\text{s}}k}} = {{\boldsymbol{T}}_{{\text{s}}k}}{{\boldsymbol{Q}}_{\text{s}}}{\boldsymbol{T}}_{{\text{s}}k}^{\text{T}}$，其中，T_mk、T_sk如式（7）.

$$ \begin{split} \left\{ \begin{array}{l} \begin{array}{l} {{\boldsymbol{T}}_{{\text{m}}k}} = \left[ {\begin{array}{*{20}{c}} {\begin{array}{*{20}{l}} {{\rm{co}}{{\rm{s}}^{\rm{2}}}(\vartheta + {T_k})}\\ {{{\sin }^{\rm{2}}}(\vartheta + {T_k})}\\ {\sin (\vartheta + {T_k})\cos (\vartheta + {T_k})} \end{array}}&{\begin{array}{*{20}{c}} {{{\sin }^{\rm{2}}}(\vartheta + {T_k})}\\ {{\rm{co}}{{\rm{s}}^{\rm{2}}}(\vartheta + {T_k})}\\ { - \sin (\vartheta + {T_k})\cos (\vartheta + {T_k})} \end{array}} \end{array}} \right. \left. {\begin{array}{*{20}{c}} { - {{\sin }^2}(\vartheta + {T_k})}\\ {{{\sin }^2}(\vartheta + {T_k})}\\ {{\rm{co}}{{\rm{s}}^{\rm{2}}}(\vartheta + {T_k}) - {{\sin }^{\rm{2}}}(\vartheta + {T_k})} \end{array}} \right] , \end{array} \\ {{\boldsymbol{T}}_{{\rm{s}}k}} = \left[ {\begin{array}{*{20}{c}} {\cos (\vartheta + {T_k})}&{\sin (\vartheta + {T_k})} \\ {\sin (\vartheta + {T_k})}&{\cos (\vartheta + {T_k})} \end{array}} \right] , \end{array} \right. \end{split} $$

(7)

式中：$ {T_k} $为曲线纤维方位角.

根据复合材料层合理论，对于n层的复合变刚度壁板本构关系可写为

$$ \begin{split} \left\{ \begin{array}{l} \left\{ {\begin{array}{*{20}{c}} {\boldsymbol{N}} \\ {\boldsymbol{M}} \end{array}} \right\} = \left[ {\begin{array}{*{20}{c}} {\boldsymbol{A}}&{\boldsymbol{B}} \\ {\boldsymbol{B}}&{\boldsymbol{D}} \end{array}} \right]\left( {\left\{ {\begin{array}{*{20}{c}} {{{\boldsymbol{\varepsilon}} _0}} \\ {\boldsymbol{\kappa}} \end{array}} \right\} + \eta \left\{ {\begin{array}{*{20}{c}} {{{\dot {\boldsymbol{\varepsilon}} }_0}} \\ {\dot {\boldsymbol{\kappa}} } \end{array}} \right\}} \right) , \\ {{\boldsymbol{F}}_{\text{s}}} = {{\boldsymbol{A}}_{\text{s}}}\left( {{\mathbf{\gamma }} + \eta {\mathbf{\dot \gamma }}} \right) , \end{array} \right. \end{split} $$

(8)

式中：${\boldsymbol{N}}、{\boldsymbol{M}}、{{\boldsymbol{F}}}_{{\text{s}}}$分别为复合变刚度壁板的膜力、弯矩、横向剪力，矩阵${\boldsymbol{A}}、{\boldsymbol{B}}、{\boldsymbol{D}}、{{\boldsymbol{A}}}_{{\rm{s}}}$分别为

$$ \begin{gathered} {\boldsymbol{A}} = \sum\limits_{k = 1}^n {{{\boldsymbol{\bar Q}}_{{\rm{m}}k}}({{\textit{z}}_k} - {{\textit{z}}_{k - 1}})} ,\;{\boldsymbol{B}} = \frac{{\text{1}}}{{\text{2}}}\sum\limits_{k = 1}^n {{{\boldsymbol{\bar Q}}_{{\rm{m}}k}}({\textit{z}}_k^2 - {\textit{z}}_{k - 1}^2)} , \\ {\boldsymbol{D}} = \dfrac{{\text{1}}}{{\text{3}}}\sum\limits_{k = 1}^n {{{\boldsymbol{\bar Q}}_{{\rm{m}}k}}({\textit{z}}_k^3 - {\textit{z}}_{k - 1}^3)} ,\;{{\boldsymbol{A}}_{\rm{s}}} = \sum\limits_{k = 1}^n {{{\boldsymbol{\bar Q}}_{{\rm{s}}k}}({{\textit{z}}_k} - {{{\textit{z}}}_{k - 1}})} . \\ \end{gathered} $$

2.3 亚音速气动力模型

气动力方面，由于高速列车马赫数较小时，气体的可压缩性可以忽略，因此，采用Dowell基于线性无黏势流理论提出的亚音速气动力模型^[15]，如式（9）.

$$ \begin{split} & \Delta p(x,t)=A(x)\left({V}_{\infty }^{2}\dfrac{{\partial }^{2}w}{\partial {x}^{2}} + 2{V}_{\infty }^{}\dfrac{{\partial }^{2}w}{\partial x\partial t} + \dfrac{{\partial }^{2}w}{\partial {t}^{2}}\right)+\\ & \quad {\left. { {D_1}(x)\left( {V_\infty ^{}\frac{{\partial w}}{{\partial t}} + V_\infty ^2\frac{{\partial w}}{{\partial x}}} \right)} \right|_{x = {{ - L/2}}}}+\\ &\quad {\left. { {D_2}(x)\left( {V_\infty ^{}\frac{{\partial w}}{{\partial t}} + V_\infty ^2\frac{{\partial w}}{{\partial x}}} \right)} \right|_{x = L/2}} \;\;\;\;, \end{split} $$

(9)

式中：$ A(x) $为壁板内气动力相互影响效应，如式（10）；$ {D}_{1}(x)、{D}_{\text{2}}(x) $为壁板边界条件对气动力的影响效应，如式（10），$ {\rho _\infty } $为空气密度.

$$ \left\{ \begin{array}{l} A(x)=-\dfrac{{\rho }_{\infty }}{{\text{π}} }{\displaystyle {\int }_{-1/2-x/L}^{1/2-x/L}\mathrm{ln}\left|\chi \right|{\rm{d}}\chi }, \\ {D}_{1}(x)= -\dfrac{{\rho }_{\infty }}{{\text{π}} }\mathrm{ln}\left(\dfrac{1}{2} + \dfrac{x}{L}\right),\\ {D}_{2}(x)=-\dfrac{{\rho }_{\infty }}{{\text{π}} }\mathrm{ln}\left(\dfrac{1}{2}-\dfrac{x}{L}\right). \end{array} \right. $$

(10)

2.4 气动弹性有限元方程

基于复合变刚度壁板黏弹性本构关系、几何关系以及气动力模型，可以写出其虚动能$ {\rm{\delta}} T $、虚应变能$ {\rm{\delta}} U $、外力虚功$ {\rm{\delta}} W $为

$$ \left\{ \begin{array}{l} \delta T = \dfrac{1}{2}\delta \displaystyle\int_{} {\rho h\left[ {{{\left( {\dfrac{{\partial u}}{{\partial t}}} \right)}^2} + {{\left( {\dfrac{{\partial \upsilon }}{{\partial t}}} \right)}^2} + {{\left( {\dfrac{{\partial w}}{{\partial t}}} \right)}^2}} \right]{\rm{d}}x{\rm{d}}y,} \\ \begin{array}{l} \delta U = \dfrac{1}{2}\delta \displaystyle\int {\left[ {\begin{array}{*{20}{c}} {{{{{\boldsymbol{\varepsilon}} }}_0}}\\ {\boldsymbol{\kappa}} \end{array}} \right]^{\rm{T}}}\left[ {\begin{array}{*{20}{c}} {\boldsymbol{A}}&{\boldsymbol{B}}\\ {\boldsymbol{B}}&{\boldsymbol{D}} \end{array}} \right]\bigg( {\left[ {\begin{array}{*{20}{c}} {{\boldsymbol{\varepsilon} _0}}\\ {\boldsymbol{\kappa}} \end{array}} \right] + } \bigg.\\ \quad\left. {\eta \left[ {\begin{array}{*{20}{c}} {{{{{\dot {\boldsymbol{\varepsilon}} }}}_0}}\\ {\dot {\boldsymbol{\kappa}} } \end{array}} \right]} \right){\rm{d}}x{\rm{d}}y + \dfrac{1}{2}\delta \displaystyle\int\limits_{} {{{\boldsymbol{\gamma }}^{\rm{T}}}{A_{\rm{s}}}\left( {{\bf{\gamma }} + \eta {\bf{\dot \gamma }}} \right){\rm{d}}x{\rm{d}}y,} \end{array} \\ \delta W = \displaystyle\int_{} {\Delta p\delta w{\rm{d}}x{\rm{d}}y} . \end{array} \right. $$

(11)

采用板单元对复合变刚度壁板进行网格划分，根据Hamilton原理，通过变分可得到壁板的有限元颤振方程为

$$ \left( {{\boldsymbol{M}} + {{\boldsymbol{M}}_{\rm{a}}}} \right){\boldsymbol{\ddot u}} + \left( {{\boldsymbol{C}} + {{\boldsymbol{C}}_{\rm{a}}}} \right){\boldsymbol{\dot u}} + \left( {{\boldsymbol{K}} + {{\boldsymbol{K}}_{\rm{a}}}} \right){\boldsymbol{u}} = {\mathbf{0}} \text{，} $$

(12)

式中：${\boldsymbol{\ddot{u}}}、{\boldsymbol{\dot{u}}}、{\boldsymbol{u}}$分别为节点加速度、速度、位移列向量；${\boldsymbol{M}}$、${\boldsymbol{C}}$、${\boldsymbol{ K }}$分别为质量矩阵、阻尼矩阵、刚度矩阵，如式（13），${{\boldsymbol{B}}}_{{\rm{b}}}、{{\boldsymbol{B}}}_{{\rm{s}}}$分别为线应变矩阵和切应变矩阵，${\boldsymbol{N}}$为形函数矩阵；$\rho $为复合材料壁板的密度；${{\boldsymbol{M}}}_{\text{a}}、{{\boldsymbol{C}}}_{\text{a}}、{{\boldsymbol{K}}}_{\text{a}}$分别为气动力产生的气动质量矩阵、气动阻尼矩阵、气动刚度矩阵，如式（14），${{\boldsymbol{N}}}_{\text{w}}$为形函数矩阵${\boldsymbol{N}}$中对应挠度自由度元素组成的矩阵.

$$ \left\{ \begin{array}{l} {\boldsymbol{M}} = \int \rho h{{\boldsymbol{N}}^{\text{T}}}{\boldsymbol{N}}{\text{d}}x{\text{d}}y, \\ {\boldsymbol{K}} = \int {{\boldsymbol{B}}_{\rm{b}}^{\text{T}}\left[ {\begin{array}{*{20}{c}} {\boldsymbol{A}}&{\boldsymbol{B}} \\ {\boldsymbol{B}}&{\boldsymbol{D}} \end{array}} \right]{{\boldsymbol{B}}_{\rm{b}}}} {\text{d}}x{\text{d}}y + \int {{\boldsymbol{B}}_{\rm{s}}^{\text{T}}{{\boldsymbol{A}}_{\rm{s}}}{{\boldsymbol{B}}_{\rm{s}}}} {\text{d}}x{\text{d}}y , \\ {\boldsymbol{C}} = \int {\eta {\boldsymbol{B}}_{\rm{b}}^{\text{T}}\left[ {\begin{array}{*{20}{c}} {\boldsymbol{A}}&{\boldsymbol{B}} \\ {\boldsymbol{B}}&{\boldsymbol{D}} \end{array}} \right]{{\boldsymbol{B}}_{\rm{b}}}} {\text{d}}x{\text{d}}y + \int {\eta {\boldsymbol{B}}_{\rm{s}}^{\text{T}}{{\boldsymbol{A}}_{\rm{s}}}{{\boldsymbol{B}}_{\rm{s}}}} {\text{d}}x{\text{d}}y , \end{array} \right. $$

(13)

$$\qquad\qquad \left\{ \begin{array}{l} {{\boldsymbol{M}}_{\rm{a}}} = \displaystyle\int {A(x)} {\boldsymbol{N}}_{\rm{w}}^{\text{T}}{{\boldsymbol{N}}_{\rm{w}}}{\text{d}}x{\text{d}}y, \\ {{\boldsymbol{C}}_{\rm{a}}} = \displaystyle\int {\left( \begin{gathered} 2A(x){V_\infty }{\boldsymbol{N}}_{\rm{w}}^{\text{T}}\frac{{\partial {{\boldsymbol{N}}_w}}}{{\partial x}} + {D_1}(x){V_\infty }{\left. {\left( {{\boldsymbol{N}}_{\rm{w}}^{\text{T}}{{\boldsymbol{N}}_{\rm{w}}}} \right)} \right|_{x = {{ - }}L/2}} + {D_2}(x){V_\infty }{\left. {\left( {{\boldsymbol{N}}_{\rm{w}}^{\text{T}}{{\boldsymbol{N}}_{\rm{w}}}} \right)} \right|_{x = L/2}} \\ \end{gathered} \right)} {\text{d}}x{\text{d}}y , \\ {{\boldsymbol{K}}_{\rm{a}}} = \displaystyle\int {\left( \begin{gathered} A(x)V_\infty ^2{\boldsymbol{N}}_{\rm{w}}^{\text{T}}\frac{{{\partial ^2}{{\boldsymbol{N}}_w}}}{{\partial {x^2}}} + {D_1}(x)V_\infty ^2{\left. {\left( {{\boldsymbol{N}}_{\rm{w}}^{\text{T}}\frac{{\partial {{\boldsymbol{N}}_{\rm{w}}}}}{{\partial x}}} \right)} \right|_{x = - L/2}} + {D_2}(x)V_\infty ^2{\left. {\left( {{\boldsymbol{N}}_{\rm{w}}^{\text{T}}\frac{{\partial {{\boldsymbol{N}}_{\rm{w}}}}}{{\partial x}}} \right)} \right|_{x = L/2}} \\ \end{gathered} \right)} {\text{d}}x{\text{d}}y . \end{array} \right. $$

(14)

3. 算法验证与结果分析

3.1 正确性验证

首先，对本文方法及程序进行正确性验证. 采用文献[12]中曲线纤维变刚度复合材料层合壁板算例，通过退化气动力相关矩阵，选取铺层方式为[0/90<0|45°>/<0|−45°>]_s的变刚度复合壁板进行了固有频率的网格收敛性分析，具体结果见表1. 分析时壁板采用四边固支约束边界条件. 为方便对比，固有频率按文献[12]进行无量纲处理.

表 1 变刚度复合壁板固有频率的网格收敛性

Table 1. Grid convergence of natural frequencies of composite panels with variable stiffness

固有频率	网格规模				文献[12]	误差/%
固有频率	5 × 5	10 × 10	20 × 20	30 × 30	文献[12]	误差/%
一阶	3.786	4.060	3.936	3.937	3.941	0.10
二阶	5.853	5.278	5.126	5.090	5.087	0.06
三阶	13.087	7.848	7.582	7.532	7.565	0.43
四阶	16.219	11.012	10.740	10.358	10.385	0.26
五阶	17.409	12.281	10.978	10.658	11.052	3.56
六阶	22.249	12.274	11.939	11.697	11.327	3.27

下载: 导出CSV

| 显示表格

从表1中可以看到：随着壁板网格规模从5 × 5增加至30 × 30，变刚度复合壁板前六阶固有频率逐渐收敛，且收敛后的各阶固有频率与文献[12]的结果最大误差不超过3.56%.

图2给出了6种铺层方式下变刚度复合壁板的固有频率，6种铺层方式分别为铺层A1[0/90<0|45°>/<0|−45°>]_s、铺层A2[0/90<15°|45°>/<−15°|−45°>]_s、铺层A3[0/90<30°|45°>/<−30°|−45°>]_s 、铺层A4[0/90<45°|40°>/<−45°|−60°>]_s、铺层A5 [0/90<45°|75° >/<−45°|−75°>]_s 、铺层A6[0/90<90°|45°>/<−45°|−90°>]_s. 壁板采用30 × 30的网格划分，四边固支约束边界条件. 从图中对比可以看出：各种铺层方式下，本文前三阶固有频率结果与文献[12]的结果吻合得都比较好. 因此，本文后续壁板稳定性分析均采用30 × 30的网格划分.

图 2 曲线纤维壁板固有频率结果对比

注：每种图例下第一个柱状条为本文解，第二个柱状条为文献[12]解.

Figure 2. Comparison of natural frequencies of curved fiber composite laminates

下载: 全尺寸图片幻灯片

其次，采用文献[16]中的两边固支铝制壁板，将本文程序退化为各向同性材料壁板进行气动弹性稳定性计算，并分析验证，结果如表2. 从表2中可以看到，本文方法计算结果与文献[16]解吻合很好.

表 2 两端固支壁板发散稳定性结果对比

Table 2. Divergence stability of the plate with bilateral fixation

厚度/ mm	板长 0.8 m		板长 0.9 m		板长 1.0 m
厚度/ mm	文献	本文	文献	本文	文献	本文
2.40	73.53	71.03	61.62	61.01	52.61	53.01
2.60	82.90	80.25	69.48	68.51	59.23	59.10
2.90	97.66	95.47	81.84	81.06	69.88	70.04
3.30	118.55	116.67	99.35	98.21	84.82	85.12

下载: 导出CSV

| 显示表格

3.2 复合变刚度壁板发散稳定性分析

1）失稳速度计算

采用复模态理论求解复合变刚度弹性、黏弹性壁板的气动弹性发散失稳速度. 图3给出了铺层方式为[0/90/<0|45>/<0|−45>]_s的弹性/黏弹性变刚度壁板在四边固支、四边简支条件下，前两阶频率随来流速度变化情况. 黏弹性材料的黏性阻尼系数η=0.002，下文中不特别注明，均取此值.

图 3 黏弹性和弹性变刚度壁板频率随气流速度变化

Figure 3. Change of the two frequencies of viscoelastic and elastic variable-stiffness panels with air velocity

下载: 全尺寸图片幻灯片

从图3中可以看到：无论变刚度壁板是弹性材料还是黏弹性材料，其一阶固有频率随着来流速度增大均是减小的，二阶固有频率随着来流速度先减小后增大；变刚度壁板材料的黏弹性效应会导致壁板二阶固有频率明显降低；当一阶固有频率减小至0时，壁板就出现发散失稳，据此可判断，四边固支的变刚度弹性、黏弹性壁板发散失稳速度分别为108 m/s和91 m/s，四边简支弹性、黏弹性变刚度壁板发散失稳速度分别为72 m/s和66 m/s. 由此可知，与弹性变刚度壁板相比，壁板材料的黏弹性效应会导致其失稳速度降低. 从图中失稳前壁板固有频率值变化剧烈程度可以看到，黏弹性壁板发散失稳出现得更突兀.

2）曲线铺设和直线铺设的比较

图4为弹性、黏弹性壁板分别采用直线和曲线纤维铺设情况下，四边固支纤维复合弹性壁板的发散稳定特性. 直线铺设与曲线铺设的铺层角考虑了7种组合方式，分别为铺层B1[0/15/45/60]、铺层B2[0/30/45/75]、铺层B3[0/15/30/60]、铺层B4[30/45/60/75]、铺层B5[45/15/30/60]、铺层B6[30/15/45/60]、铺层B7[45/30/60/75]. 其中，曲线纤维铺设时，纤维方向角T₀=0°，T₁=45°. 从图中可以看到：无论是弹性壁板还是黏弹性壁板，在纤维曲线铺设下的失稳速度相比直线纤维铺设都有提高. 采用铺层B7时，曲线纤维壁板失稳速度相比直线纤维壁板可提高将近50%.

图 4 直线铺设与曲线铺设对变刚度壁板发散临界速度的影响

Figure 4. Comparison of variable-stiffness panel instability with straight and curved fiber laying

下载: 全尺寸图片幻灯片

3）铺层主方向角的影响

图5为铺层主方向对曲线纤维变刚度壁板发散临界速度的影响. 7种不同的铺层方式为铺层C1[0/0]_s、铺层C2[0/15]_s、铺层C3[0/30]_s、铺层C4[0/45]_s、铺层C5[0/60]_s、铺层C6[0/75]_s、铺层C7[0/90]_s等，壁板采用四边固支约束. 从图中可以看出，随着铺层主方向角从0°~90° 增大，无论是变刚度弹性壁板还是变刚度黏弹性壁板，壁板发散失稳速度呈增大的趋势. 由此可见，通过增大铺层主方向角可以提高复合材料变刚度壁板发散速度.

图 5 铺层主方向对曲线纤维变刚度壁板发散临界速度影响

Figure 5. Influence of ply direction on divergence critical speed of variable-stiffness panels with curved fibers

下载: 全尺寸图片幻灯片

4）曲线纤维方向角的影响

进一步研究曲线纤维方向角对弹性、黏弹性变刚度复合壁板发散稳定性的影响. 壁板为四边固支约束，其铺层形式为[0/90/<T₀|T₁>/<−T₀|−T₁>]_s. 考虑T₀和T₁单独变化的两种情况. 图6为T₁不变，T₀=−90°~90°（间隔15°）变刚度壁板发散失稳特性随T₀的变化规律. 图7为T₀不变，T₁= −90°~90°（间隔15°）变刚度壁板发散失稳特性随T₁的变化规律.

图 6 曲线纤维变刚度壁板失稳速度随T₀的变化（T₁不变）

Figure 6. Changes of divergence characteristics with T₀ for variable-stiffness plate with curved fibers (T₁ unchanged)

下载: 全尺寸图片幻灯片

图 7 曲线纤维变刚度壁板失稳速度随T₁的变化（T₀不变）

Figure 7. Changes of flutter characteristics with T₁ for variable-stiffness plate with curved fibers (T₀ unchanged)

下载: 全尺寸图片幻灯片

从图6中可以看出：当T₁不变、T₀变化时，变刚度弹性、黏弹性壁板发散速度随着T₀的增大均呈先减小后增大的趋势，呈“凹”形分布；对变刚度弹性壁板，在T₀从−90°~−60° 时，壁板在亚音速已不发生失稳，T₀从15°~45° 时，壁板发散失稳速度出现极小值范围，基本关于T₀ = 30° 对称分布；对于变刚度黏弹性壁板，其发散失稳速度分布规律与弹性壁板类似，只是变刚度黏弹性壁板发散失稳速度极小值相比变刚度弹性壁板的小.

从图7中可以看出：当T₀不变、T₁变化时，变刚度弹性、黏弹性壁板壁板发散速度随着T₁从−90°～90° 也呈“凹”形分布；壁板发散失稳速度关于T₀=0° 对称分布，发散失稳速度变化相对平缓.

由此可见，复合变刚度壁板曲线纤维参数T₀、T₁对变刚度弹性、黏弹性壁板发散稳定性均具有明显影响，而且T₁影响更为明显一些. 此外，曲线纤维复合变刚度壁板在发散稳定性方面存在较大的优化空间，通过调整曲线纤维的路径可以提高复合材料变刚度壁板的发散速度.

5）黏弹性材料阻尼系数的影响

图8给出了黏弹性材料阻尼系数对曲线纤维复合壁板发散稳定性的影响. 从图可以看到：随着黏弹性材料阻尼系数增大，曲线纤维复合壁板的失稳速度减小.

图 8 黏弹性阻尼对曲线纤维变刚度壁板发散临界速度影响

Figure 8. Influence of viscoelastic damp on divergence critical speed of variable-stiffness panels with curved fibers

下载: 全尺寸图片幻灯片

6）边界条件影响

图9给出了不同约束情况下，变刚度弹性、黏弹性壁板的发散稳定性情况. 壁板铺层形式为[0/90/<T₀|T₁>/<−T₀|−T₁>]_s. 从图中可以看出：约束条件对变刚度壁板失稳速度影响很明显，边界约束越强，壁板的稳定性越好，发散失稳临界速度越大.

图 9 约束对曲线纤维变刚度壁板发散失稳临界速度影响

Figure 9. Influence of boundary conditions on divergence critical speed of variable-stiffness panels with curved fibers

下载: 全尺寸图片幻灯片

4. 结　论

本文建立了亚音速流场中曲线纤维变刚度弹性、黏弹性壁板气动弹性稳定性分析模型，与文献[12]、文献[16]结果进行对比，验证了本文模型及算法的正确性，进而采用复模态理论求解复合变刚度壁板发散临界速度，在此基础上进行了关键参数影响分析，得出的主要结论如下：

1）变刚度壁板材料的黏弹性效应会导致其失稳速度降低，与弹性壁板相比，黏弹性壁板失稳前固有频率值变化更剧烈，发散失稳出现得更突兀；

2）通过调整铺层主方向角可以实现复合材料变刚度弹性、黏弹性壁板的发散速度的提高，通过调整曲线纤维的路径也可以提高复合材料变刚度壁板的发散速度，因此，曲线纤维复合变刚度壁板在发散稳定性方面存在较大的优化设计空间；

3）弹性、黏弹性壁板发散速度随着曲线纤维路径参数T₀、T₁从 −90° ~ 90° 变化呈现“凹”形分布，且参数T₁对壁板发散稳定性的影响相比T₀更明显；

4）边界条件对曲线纤维复合材料弹性、黏弹性壁板的发散临界速度有较大影响，边界约束越强，壁板的稳定性越好，发散失稳临界速度越大.

致谢：民航航空器适航审定技术重点实验室开放基金（SH2020112705）.

图 1 曲线纤维层合板

Figure 1. Curvilinear fibers composite panels

下载: 全尺寸图片幻灯片

图 2 曲线纤维壁板固有频率结果对比

注：每种图例下第一个柱状条为本文解，第二个柱状条为文献[12]解.

Figure 2. Comparison of natural frequencies of curved fiber composite laminates

下载: 全尺寸图片幻灯片

图 3 黏弹性和弹性变刚度壁板频率随气流速度变化

Figure 3. Change of the two frequencies of viscoelastic and elastic variable-stiffness panels with air velocity

下载: 全尺寸图片幻灯片

图 4 直线铺设与曲线铺设对变刚度壁板发散临界速度的影响

Figure 4. Comparison of variable-stiffness panel instability with straight and curved fiber laying

下载: 全尺寸图片幻灯片

图 5 铺层主方向对曲线纤维变刚度壁板发散临界速度影响

Figure 5. Influence of ply direction on divergence critical speed of variable-stiffness panels with curved fibers

下载: 全尺寸图片幻灯片

图 6 曲线纤维变刚度壁板失稳速度随T₀的变化（T₁不变）

Figure 6. Changes of divergence characteristics with T₀ for variable-stiffness plate with curved fibers (T₁ unchanged)

下载: 全尺寸图片幻灯片

图 7 曲线纤维变刚度壁板失稳速度随T₁的变化（T₀不变）

Figure 7. Changes of flutter characteristics with T₁ for variable-stiffness plate with curved fibers (T₀ unchanged)

下载: 全尺寸图片幻灯片

图 8 黏弹性阻尼对曲线纤维变刚度壁板发散临界速度影响

Figure 8. Influence of viscoelastic damp on divergence critical speed of variable-stiffness panels with curved fibers

下载: 全尺寸图片幻灯片

图 9 约束对曲线纤维变刚度壁板发散失稳临界速度影响

Figure 9. Influence of boundary conditions on divergence critical speed of variable-stiffness panels with curved fibers

下载: 全尺寸图片幻灯片

表 1 变刚度复合壁板固有频率的网格收敛性

Table 1. Grid convergence of natural frequencies of composite panels with variable stiffness

固有频率	网格规模				文献[12]	误差/%
固有频率	5 × 5	10 × 10	20 × 20	30 × 30	文献[12]	误差/%
一阶	3.786	4.060	3.936	3.937	3.941	0.10
二阶	5.853	5.278	5.126	5.090	5.087	0.06
三阶	13.087	7.848	7.582	7.532	7.565	0.43
四阶	16.219	11.012	10.740	10.358	10.385	0.26
五阶	17.409	12.281	10.978	10.658	11.052	3.56
六阶	22.249	12.274	11.939	11.697	11.327	3.27

下载: 导出CSV

表 2 两端固支壁板发散稳定性结果对比

Table 2. Divergence stability of the plate with bilateral fixation

厚度/ mm	板长 0.8 m		板长 0.9 m		板长 1.0 m
厚度/ mm	文献	本文	文献	本文	文献	本文
2.40	73.53	71.03	61.62	61.01	52.61	53.01
2.60	82.90	80.25	69.48	68.51	59.23	59.10
2.90	97.66	95.47	81.84	81.06	69.88	70.04
3.30	118.55	116.67	99.35	98.21	84.82	85.12

下载: 导出CSV

参考文献(16)

[1]	GŰRDAL Z, TATTING B F, WU C K. Variable-stiffness composite panels:effects of stiffness variation on the in-plane and buckling response[J]. Composites:Part A:Applied Science and Manufacturing, 2008, 39(5): 911-922. doi: 10.1016/j.compositesa.2007.11.015
[2]	HYPER M W, CHARETTE R F. Use of curvilinear fiber format in composite structure design[J]. AIAA Journal, 1991, 29(6): 1011-1015. doi: 10.2514/3.10697
[3]	LOPES C S, GÜRDAL Z, CAMANHO P P. Variable-stiffness composite panels:Buckling and first-ply failure improvements over straight-fibre laminates[J]. Computers & Structures, 2008, 86(9): 897-907. doi: 10.1016/j.compstruc.2007.04.016
[4]	HAMED A, PEDRO R, DEMOURA M F S F. Large deflection and stresses in variable stiffness composite laminates with curvilinear fibres[J]. International Journal of Mechanical Sciences, 2013, 73: 14-26. doi: 10.1016/j.ijmecsci.2013.03.013
[5]	马洪涛. 变刚度复合材料层合板的力学性能研究[D]. 哈尔滨: 哈尔滨工业大学, 2012.
[6]	聂国隽,朱佳瑜. 纤维曲线铺放的复合材料层合板的自由振动分析[J]. 力学季刊,2016,37(2): 274-283. NIE Guojun, ZHU Jiayu. Free vibration analysis of composite laminates with curvilinear fibers[J]. Chinese Quarterly of Mechanics, 2016, 37(2): 274-283.
[7]	马成. 复合材料纤维曲线铺放层合板减振性能分析[D]. 南京: 南京航空航天大学, 2019.
[8]	GROH R M J, WEAVER P M. Buckling analysis of variable angle tow,variable thickness panels with transverse shear effects[J]. Composite Structures, 2014, 107: 482-493. doi: 10.1016/j.compstruct.2013.08.025
[9]	HAO P, YUAN X J, LIU C, et al. An integrated framework of exact modeling,isogeometric analysis and optimization for variable-stiffness composite panels[J]. Computer Methods in Applied Mechanics and Engineering, 2018, 339: 205-238. doi: 10.1016/j.cma.2018.04.046
[10]	孙士平,张冰,邓同强,等. 复合载荷作用变刚度复合材料回转壳屈曲优化[J]. 复合材料学报,2019,36(4): 1052-1061. SUN Shiping, ZHANG Bing, DENG Tongqiang, et al. Buckling optimization of variable stiffness composite rotary shell under combined loads[J]. Acta Materiae Compositae Sinica, 2019, 36(4): 1052-1061.
[11]	VAHID K, JAMSHID F. Supersonic panel flutter of variable stiffness composite laminated skew panels subjected to yawed flow by using NURBS-based isogeometric approach[J]. Journal of Fluids and Structures, 2018, 82: 198-214. doi: 10.1016/j.jfluidstructs.2018.07.002
[12]	欧阳小穗,刘毅. 高速流场中变刚度复合材料层合板颤振分析[J]. 航空学报,2018,39(3): 221539. OUYANG Xiaosui, LIU Yi. Panel flutter of variable stiffness composite laminates in supersonic flow[J]. Acta Aeronautica et Astronautica Sinica, 2018, 39(3): 221539.
[13]	TOURAJ F, DAVOOD A, HASAN K. Flutter improvement of a thin walled wing-engine system by applying curvilinear fiber path[J]. Aerospace Science and Technology, 2019, 93(2): 105353.1-105353.14.
[14]	许家宝. 亚麻/玻璃纤维混杂复合材料的动态粘弹性及湿热性能研究[D]. 南京: 南京航空航天大学, 2019.
[15]	KORNECKI A, DOWELL E H, O’BRIEN J. On the aeroelastic instability of two-dimensional panels in uniform incompressible flow[J]. Journal of Sound and Vibration, 1976, 47(2): 163-178.
[16]	GUO Y, LI F M. Chaotic motion of a composite laminated plate with geometric nonlinearity in subsonic flow[J]. International Journal of Non-Linear Mechanics, 2013, 50: 81-90. doi: 10.1016/j.ijnonlinmec.2012.11.010