热环境下复合材料壁板的非线性颤振特性分析

摘 要：复合材料因其优异的性能已被广泛应用于飞行器的结构中。在高超声速飞行时，由于非定常气动载荷和热效应的共同作用，飞行器复合材料壁板会出现屈曲及非线性颤振现象。本文基于均匀温度场建立了考虑热效应的复合材料薄壁结构的非线性有限元模型，结合气动-结构二阶松耦合算法，针对复合材料壁板在Ma∞=5时的非线性气动弹性响应进行分析，探讨了复合材料铺层方向、温度载荷等对壁板颤振特性的影响。结果表明，壁板响应受到温度载荷、铺层方向及来流动压等多方面的影响。斜交壁板的颤振动压高于正交壁板颤振动压（增量大于55%）;随着温度载荷的增大，壁板颤振动压会下降 ，最大降幅约为动压的80%。当出现热屈曲时，颤振动压又会缓慢增大;当温度超过结构的临界屈曲温度时，壁板在亚临界条件下会出现多个屈曲位置，随着动压的变化而发生转变。

关键词：高超声速;壁板颤振;复合材料;热效应;非线性气动弹性;热屈曲;飞行器

中图分类号：TJ760; V214.8

文献标识码：A

文章编号：1673-5048（2022）04-0091-09

DOI：10.12132/ISSN.1673-5048.2021.0183

0 引 言

为了减轻结构重量，改善结构性能，各种新型复合材料被大量应用于高超声速飞行器结构中^[1-2]。飞行器高超声速飞行时，会出现气动加热，复合材料薄壁结构在气动载荷和热载荷作用下发生变形，导致流场的边界发生改变，使得作用在结构上的载荷发生变化，进一步影响薄壁结构的响应，使其产生新的变形或者振荡^[3]。这种结构、气动和热之间的相互耦合作用，造成了壁板的热气动弹性响应问题^[4]。与传统翼、舵等升力面构型不同的是，薄壁结构在非定常载荷作用下会表现出非线性流-固耦合特征：其位移响应通常与厚度为相同量级，面内应力存在，弯曲和拉伸之间会发生耦合，在周围受约束条件下，形变一般不会引发结构的迅速破坏，表现为几何非线性^[5-7]。

国内外一些学者在该领域有一些很有意义的研究，Nydick等^[8]基于Marguerre薄壳理论和Galerkin方法求解结构变形方程，分析了高超声速气流中三维正交各向异性薄板的极限环幅值和振荡形态，指出壁板的气动弹性响应受到初始曲率、气动热、激波运动等的影响。Shiau等^[9]研究了高超声速气流中层合板的非线性颤振响应，指出材料各向异性对壁板的极限环和混沌运动形态等有显著影响。Gray^[10]结合Von-Karman板理论和三阶活塞理论分析了复合材料薄壁板的非线性气动弹性响应，分析了支承条件、系统参数对颤振特性的影响。Kouchakzadeh等^[11]基于Galerkin方法计算了超声速气流中复合材料层合板的颤振特性，并讨论了面内应力、纤维方向、气动阻尼等对非线性响应的影响。

梁路等^[12]研究了具有复合材料壁板结构的机翼气动弹性问题，结果表明复合材料不同铺层的单层厚度比对机翼的刚度影响很大。欧阳小穂等^[13]研究了高速气流下变刚度复合材料壁板的非线性颤振响应，分析了支承条件和纤维方向等对颤振边界的影响。林华刚^[14]研究了非均匀温度场中复合材料圆柱壳的颤振特性，并进行了气动热-气动弹性耦合的计算分析。苑凯华等^[15]分析了复合材料壁板的非线性颤振特性，结果表明高超声速气流中气动力的非线性和热载荷对壁板颤振的振动幅值的影响更为明显。杨智春等^[16-17]利用有限元方法分析了不同铺层方向和铺层顺序对层合复合材料壁板颤振特性的影响，发现随温度升高，发生频率重合型颤振的耦合模态会发生演变，使得壁板有可能在较低的速压下发生颤振。

可以看出，在复合材料壁板的非线性气动弹性的研究中，大都采用了Galerkin方法和活塞理论。需要指出的是，活塞理论求解非定常气动载荷的适用性有一定限制（2<Ma∞<5），且无法处理流场的三维效应和黏性效应。而基于Galerkin方法的结构变形分析中，其只适合于驻波振荡及特定支撑条件。

对于大尺度高超声速飞行器而言，在持续稳定飞行状态下，飞行器的表面存在近似均匀升温区域。本文基于均匀温度场针对复合材料薄壁结构在高超声速气流中的非线性气动弹性响应展开分析，建立了考虑热效应的复合材料壁板的非线性有限元求解模型，并将其应用于气动-结构的二阶松耦合方法中，用于分析考虑非线性非定常气动载荷和温度载荷共同作用下的非线性气动弹性问题，揭示非线性耦合的基本关系。航空兵器 2022年第29卷第4期

屈佑文，等：热环境下复合材料壁板的非线性颤振特性分析

1 复合材料薄壁结构的非线性结构动力学建模

1.1 复合材料薄壁壳元的切线刚度矩阵

薄壁结构中面变形可以看作是薄膜变形和弯曲变形相互耦合的结果，引入一阶剪切变形假设。复合材料薄壁结构如图1所示，其中（X， Y， Z）为单元参考坐标系，（L， T）为铺层材料参考坐标系。

考虑温度场的变化[100 ℃，900 ℃]，基于本文的热效应有限元模型进行模态分析，并与Nastran的结果进行比较。图5给出了在不同温度下前6阶热模态频率的计算结果对比。可以看出，本文的计算结果与Nastran的结果颇为接近，除了第4阶模态（该模态为扭转模态）有稍许差异外，其他模态的频率吻合较好。第4阶模态最大误差不超过9%，这种差异可能与两种算法的非线性刚度矩阵计算的不同有关。

由图可知，随温度的升高，各阶频率会随之降低，更高阶的频率下降幅度更大。这种特性很容易造成各阶模态运动之间的耦合，从而诱发如壁板颤振等不利的气动弹性现象。前两阶频率随温度出现了先降低后增高的趋势。以第1阶结构模态为例，其频率随温度的升高先变小，在400 ℃时达到最小，在超过429.5 ℃时，1阶频率又上升，此时结构发生了热屈曲，热屈曲以后会产生较大的变形。屈曲对板频率的影响，其原因较为复杂，与薄板屈曲构型、温度梯度造成的热应力以及振动模态有关。本文计算中，基频模态在屈曲后频率随温度升高而增大。这也与Schneider等^[23]的研究结果相符。

4.3 非定常气动力计算验证

以AGARD CT5标模为例进行验证，翼型为NACA0012，马赫数为0.755，翼型非定常运动规律为α（t）=0.016+2.51·sin（2×0.081 4 t） 。图6显示了本文计算的升力系数和俯仰力矩系数与实验值的比较。两者吻合很好，说明了本文所述方法应用于无黏非定常流场计算时的可靠性。

4.4 三维壁板颤振分析验证

考虑一个三维各向同性壁板，如图7所示。其几何和材料特性为：a/b=1，h/a=0.002，μ=0.3，μs=ρaa/（ρsh）=0.1，颤振分析时选取的计算状态为：超声速Ma∞=1.2， 壁板的流场网格为H型，包括121×121×39节点，结构有限元模型由1 152个三角壳元构成。

图8显示了壁板极限环振荡幅值和频率，并与Dowell^[6]和Gordnier^[24]的结果进行了对比。显然，与Gordnier^[24]的结果相符较好。Gordnier的求解器中同样采用了基于CFD的非定常载荷计算，而Dowell^[6]幅值稍高，因为其选用了线性势流理论来计算流场。这个算例也说明了本文耦合计算方法在处理气动力非线性和几何大变形结构非线性引起的气动弹性响应上的有效性。

5 壁板的非线性热气动弹性响应分析

5.1 壁板模型与参数

壁板几何和材料属性参数为：a=b=1.0 m，h=0.005 m，EL/ET=40.0，GLT/ET=0.6，GTT/ET=0.5，μLT=0.25，ρs=1 500，质量比μs=ρaa/（ρsh）=0.2。壁板四个边界均固支，壁板包含5个等厚度铺层。考虑两种铺层：正交铺层板[0°/90°/0°/90°/0°]（cross-ply）和斜交铺层板[45°/-45°/45°/-45°/45°]（angle-ply）。计算工况为Ma∞=5，给定温度条件分别为ΔT=0 ℃， 50 ℃， 100 ℃， 200 ℃，分析壁板在不同动压下的非线性热气动弹性响应，以x/a=0.75， y/b=0.5作为参考点。

5.2 网格无关性和时间步长收敛性分析

为了分析本文数值算法中流场网格对气动弹性响应的无关性，以及选取的耦合计算时间步长的收敛性影响，针对正交铺层在Ma∞=5， λ=500（λ=ρ∞V²∞a³/D，其中D为无量纲刚度）， ΔT=50 ℃展开分析。

壁板流场网格为H型，考核了三种网格密度，沿壁板顺气流方向、沿展向以及法向的网格点分别设置为：81×81×29（壁板表面41×41）、121×121×39（壁板表面61×61）和181×181×49（壁板表面81×81）。选取三个不同的无量纲时间步长：Δτ=0.015， 0.03， 0.06。采用Richardson外插方法^[25]考核网格和时间步长的影响，图9显示了以归一化无量纲振荡幅值为考核的收敛情况。

综合考虑计算精度和效率，本文后续计算的网格选取为121×121×39，耦合计算时间步长Δτ=0.03。图10显示了壁板流场网格，壁板在下表面中心位置，下表面设置为无穿透条件，流入边界上的变量由自由流超声速条件确定。在流出边界上，通过外推获得变量，其他边界变量由远场黎曼不变边界条件定义。图11显示了y/b=0.5剖面处壁板发生变形后的网格图。壁板流-固耦合响

应计算中，不同λ下的流场初始值取的是刚性壁板的定常流场。

5.3 正交铺层壁板在不同温度下的响应分析

（1） ΔT = 0 ℃下的响应分析

对于正交壁板，ΔT=0 ℃时，发现λ=766时出现了颤振，其相平面图如图12（a）所示。随着动压的增大，后颤振的幅值和频率也会随之增加，在较大动压下，相图会出现粉刺现象，与高阶模态介入有关， 如图12（b）所示。图13～14分别显示了λ=1 500、参考点振荡达到ф =0°（正向极值）和ф =180°（负向极值）时的壁板变形云图和流场压强云图。可以看出，此时的振荡形式表现除了结构的2阶模态之外，还有高阶形态参与其中。从压强云图来看，壁板参考点正向极值时前缘有激波，后缘为膨胀波;其负向极值时前缘有膨胀波，后缘有激波出现。

λ=1 500 at ΔT=0 ℃

（2） ΔT=50 ℃下的响应分析

当ΔT=50 ℃时，在λ=377时出现了颤振，颤振动压相比未施加温度载荷的结果变小约50%，其相平面图如图15（a）所示。随动压增大，颤振幅值和频率也会随之增加， 如图15（b）所示。

（3）ΔT=100 ℃下的响应分析

当ΔT=100 ℃时，壁板在亚临界下会出现静态屈曲。在不同动压条件下，该屈曲位置由正向朝负向发生了转变， 如图16所示。当动压增大到λ=120时出现了低频颤振，此时壁板无法保持在热屈曲的平衡位置，受