叶柳青, 叶正寅, 洪 正, 叶 坤

(西北工业大学 航空学院，西安 710072)

吸气式超声速/高超声速飞行器因具有较高的潜在军事与民用价值，已经成为了各国在航空航天领域投入巨大人力物力来研发的热点和关键性技术之一[1-3]。而冲压/超燃冲压发动机[4-7]作为此类飞行器推进系统的核心，在一定程度上决定了飞行器最终能达到的技术高度。冲压/超燃冲压发动机是一类结构简单的吸气式发动机，它们直接利用空气作为燃料氧化剂，使得飞行器具有连续可控的推力，具备更高的飞行灵活性和经济性。在冲压/超燃冲压发动机内部存在激波、膨胀波等复杂的波系[8-10]，而满足结构质量和主动冷却的设计要求需要采用薄壁结构, 因此激波主导流动中弹性壁板的动力学问题是吸气式高速飞行器结构设计与优化的重要基础科学问题。

近二十年来，许多学者针对激波主导流动中壁板的气动弹性问题展开研究，取得了一系列研究成果。通过Von-Karman大变形板弯理论考虑壁板几何非线性，并且根据CFD(computational fluid dynamics)代理模型和热传导理论，Miller等[11]建立了基于流固热的气动弹性分析模型，该模型可以对激波作用下的壁板气动弹性响应进行长时间的计算。假定激波位置固定在二维壁板的中点不变，Visbal[12-13]通过Euler和Navier-Stokes方程来分别求解变形壁板上的无黏和黏性流场，利用双向流固耦合的数值方法研究了二维弹性壁板在激波作用下的动力学特性，结果表明当激波强度较小时，气动弹性稳定性显著提高；而当激波强度较大时，临界颤振动压远小于没有激波作用时壁板的临界颤振动压，稳定性明显降低。Boyer等[14]将Visbal的二维气弹模型拓展到三维，研究了三维无黏流场中在激波的作用下弹性壁板的动力学响应，结果表明强激波的作用会增大壁板振动响应的幅值和频率，而弱激波的作用将显著增大临界颤振动压，即提高壁板的气动弹性稳定性。

Brouwer[15]提出了可用当地活塞流理论来预测激波作用下弹性壁板表面的非定常气动力。基于Von-Karman大变形板弯理论和当地活塞流理论，Ye等[16-17]建立了激波作用在二维壁板的任意位置时受热壁板的气动弹性稳定的理论分析模型，推导出了壁板发生气动弹性失稳的边界条件，并分析了激波强度、激波冲击位置对临界颤振动压、振动响应幅值和频率的影响。结果表明，随着激波冲击位置从壁板一端向另一端移动变化时，临界颤振动压、振动幅值和频率都不是单调变化的。采用双向流固耦合算法，李映坤等研究了无黏流场中斜激波冲击作用下二维曲壁板的气动弹性响应特性，着重分析了弯曲高度和动压对系统振动特性的影响。结果显示较小的弯曲高度可降低壁板颤振临界动压，而当弯曲高度进一步增大后，临界颤振动压迅速提高，并且准周期无规则运动状态被激发出来。但是，这些研究都假定激波位置固定不变，目前仅有少量的文献考虑激波位置随时间发生改变。Dennis[18]采用风洞试验对在快速移动的激波作用下弹性壁板的气动弹性特性进行了初步的研究，结果表明，与激波位置固定不变的情况相比，弹性壁板在快速移动的激波作用下振动幅值显著提高。

国内外学者对板的非线性共振问题已开展了较多的研究工作。在20世纪90年代，Okajima等[19]对悬臂矩形板的共振频率与曲率和厚度变化的关系进行了研究。在均匀分布的载荷作用下，Du等[20]对阻尼夹层圆板的非线性超谐共振问题进行了研究。基于Von-Karman 大变形理论，Xue等[21]建立了磁电弹性薄板的非线性无阻尼强迫振动的数学模型，并采用改进的L-P(lindstedt-poincare)法对问题进行求解，详细地分析了板的厚度、外激励力、压电材料、压磁材料以及磁电弹材料等对系统主共振区间、弹簧刚度特性和振幅跳跃现象的影响。基于多尺度法，张小广等[22-24]对矩形板和圆板的主共振问题进行了大量的研究。在矩形板方面，张小广等对四边固支约束的陶瓷-金属材料功能梯度矩形板的主共振问题进行了理论研究。在圆板方面，Hu等研究了静磁场作用下导电旋转圆板的主共振问题，并详细分析了磁感性强度、转速以及静载等参数对频率和振幅的影响。近年来，马冰冰等[24]对在常磁场引起的静载荷作用下铁磁圆板的主共振问题进行了研究，并分析了不同电磁参量对共振振幅的影响。采用一种简化的气动模型[25]对作用于板上的气动载荷进行描述，李文强等[26]对气动载荷作用下旋转运动导电圆板的主共振问题进行了研究。以上关于壁板主共振问题的研究工作中并没有考虑激波的作用。

实际上，冲压/超燃发动机内部的激波串往往存在振荡特征[27-29],这种激波串振荡会给发动机壁面带来严酷的压力脉动载荷。就作者目前所知，比较缺乏对振荡激波作用下壁板的主共振特性的研究。本文基于Von-Karman 大变形理论和当地活塞流理论，采用 Galerkin 方法建立了振荡激波作用下壁板振动的动力学模型，通过龙格-库塔法对非线性动力学方程进行数值计算获得系统非线性振动响应，详细分析了振荡激波强度、激波振荡幅值、温度、激波振荡的中心位置、来流马赫数对系统幅频响应特性的影响，尤其是振动突跳与双稳态特性。

1 振荡激波作用下受热壁板的非线性动力学方程

1.1 控制方程

当三维平壁板在流向上的尺寸远小于另一个方向的尺寸时，它可以简化为二维平壁板模型，这里仅考虑二维各向同性材料的壁板。振荡激波流场中的二维平壁板，边界条件为两端简支，如图1所示。壁板的长度为l，厚度为h，其中壁板厚度远小于壁板长度。壁板的上表面受振荡激波流场中非定常气动载荷的作用，入射斜激波前的气流密度、速度和马赫数分别为ρu,l，Uu,l，Mau,l，反射斜激波后的气流密度、速度和马赫数分别是ρu,r，Uu,r，Mau,r。壁板下表面受空腔压力的作用，空腔压力为pd。基于Hamilton原理和Von-Karman大变形板弯理论，壁板在振荡激波作用下的气动弹性方程[30]为

(1)

1.2 热应力

在飞行器高马赫数(一般大于2.2)飞行时，气动加热的影响一般不容忽视，还应该考虑由于气动加热产生的热应力[31-33]。假设壁板受热达到稳态并且温度均匀分布，设T0为初始温度，由温升ΔT=T-T0引起的面内热应力为

(2)

式中，α为材料的热膨胀系数。

1.3 气动力

在分析振荡激波作用下受热壁板的主共振特性中，能够有效地预测激波振荡流场中的非定常气动载荷是至关重要的。活塞理论在20世纪50年代被Lighthill[34]和Ashely[35]等提出后，由于具有简洁的表达式，且具有较好的实用性，被广泛应用于超声速、高超声速非定常气动力的计算。该理论假设当马赫数(Ma)≫1时，机翼所产生的扰动沿翼面法向传播，而机翼表面各点间的相互影响很小，就像气缸中活塞所产生的的扰动传播一样。在等熵假设条件下，通过动量定理和等熵关系式可以得到经典活塞理论计算公式为

(3)

式中：P为翼面表面任意一点的压力；P∞为来流压强；a∞为来流音速；γ为空气比热比；νn为翼面法向速度，具体表达式如下

(4)

式中，U∞为来流速度。

当来流马赫数较小时，一般采用修正的线性活塞理论[36]

(5)

Brouwer等已经证实了当地活塞流理论预测激波主导流场中气动力的可行性。在当地活塞流理论中，相应的自由来流流动参数用当地参数来代替。本文研究振荡激波作用下弹性壁板的主共振特性，主要考虑结构非线性，因此这里采用当地1阶活塞流理论，得到壁板上表面的气动力如下。

入射斜激波前的气动载荷为

(6)

反射斜激波后的气动载荷为

(7)

(8)

反射斜激波后壁板的上下压差为

(9)

1.4 无量纲化

引入无量纲参数

将上述无量纲参数代入式(1)，则无量纲化后的控制方程为

(10)

由于壁板两端简支，相应的无量纲边界条件为

(11)

实际上，本文将壁板两端的约束简化为固定铰支座，如图1所示。因此除了式(11)中所表现的边界条件(即壁板的两端可以转动，并且沿着垂直方向不可以移动)以外，壁板两端沿着水平方向也是不可以移动的。

图1 受振荡激波作用的二维壁板

1.5 Galerkin离散

采用Galerkin方法对式(10)进行离散，将位移函数展开成各阶谐波模态的叠加。对于简支边界条件，位移函数为

(12)

将式(12)代入到式(10)，对方程各项同时乘以sin(jπx)(j=1,2,…,∞)，并沿着板长积分。

这里假设壁板上表面的激波将以壁板上某一点为中心作简谐振动，入射斜激波冲击在壁板上的位置记为ξk(ξk为无量纲量，ξk=xk/l)，则ξk=ξ0+Acos(ωτ)，ξ0为中心位置，A为激波振荡的幅值，ω为激波振荡的频率。由于振荡激波的存在，壁板左右两端压差不同，采用Galerkin方法对气动力项处理如下

(13)

(14)

其中

Q=Q1+Q2+Q3+Q4+Q5，

2 模型验证

关于受热壁板在振荡激波作用下的主共振问题，在第1章中已经建立了相应的动力学分析模型。本文采用4阶Runge-Kutta法对此动力学模型(即式(14))进行数值积分，取无量纲时间步长为Δτ=0.001，初值为W0(ξ)=1，观察壁板上顺气流75%长度(ξ=0.75)处的时间历程响应，分析系统的幅频响应特性，从而探讨振荡激波作用下受热壁板的主共振特性。

表1 不同激波角下的静压比和动压比

2.1 数值方法验证

图2 本文结果与Dowell的研究结果比较

图3 本文结果与Ye等的研究结果比较

2.2 模态收敛性分析

为了选择合适的模态数目求解振荡激波作用下受热壁板的非线性气动弹性响应，有必要研究不同模态数目对响应结果的影响，即：收敛性分析。选取激波角为σ=18°，假定激波围绕壁板中点作简谐振动，即ξ0=0.500，激波振荡幅值为A=0.100，振荡频率为ω=10，图4显示模态数目分别取N=2、4、6时，振动幅值随动压的变化。其中：WA为弹性壁板的振动幅值。从图4可知，用2阶谐波模型(N=2)算得的幅值与用4、6阶谐波模型(N=4、6)算得的幅值基本一致。由此可见，当来流动压较小时，受热壁板在振荡激波作用下作强迫振动，描述这种情况下的壁板变形至少取2阶谐波模态。为保证计算的精度，在下面所有的分析中，均采用6阶谐波模型来进行计算。

图4 模态数目对振动幅值的影响

3 振荡激波作用下受热壁板主共振特性分析

采用4阶Runge-Kutta法对系统动力学方程进行数值计算，得到系统的幅频响应特性曲线，探讨振荡激波作用下受热壁板的主共振特性，尤其是振动突跳与双稳态特性。其中，幅频响应特性曲线分别采用升频扫描和降频扫描计算[37]获得，即采用前一个外激励频率激励下响应的稳态解作为下一个外激励频率计算的初始条件进行数值计算。本文重点关注主共振区，因此主要截取壁板1阶固有频率附近的幅频响应曲线进行分析。壁板的第j阶无量纲固有频率[38]为ωj=(jπ)2。影响振荡激波作用下受热壁板主共振特性的因素有很多，这里分别详细地分析了振荡激波强度、激波振荡幅值、温度、激波振荡的中心位置、来流马赫数以及下壁面背压对壁板主共振特性的影响规律。

3.1 振荡激波强度对受热壁板主共振特性的影响

首先分析振荡激波强度的影响，固定来流马赫数为Ma=3.5，来流动压为λu,l=30，选取不同的激波角σ为17°、18°、20°、22°、24°，由斜激波关系式计算得到的入射斜激波前与反射斜激波后的无量纲动压比和无量纲静压比见表1。假定激波围绕着壁板中点作简谐振荡，即ξ0=0.500，激波振荡的幅值为A=0.010。温升为ΔT/ΔTcr=0.1。采用4阶Runge-Kutta法求解系统动力学方程，得到不同振荡激波强度下受热壁板的幅频响应特性曲线，如图5所示。

其中，图5实线为激波振荡频率不断增大的正向扫频曲线，虚线为激波振荡频率不断减小的反向扫频曲线，后文幅频响应曲线中实线和虚线的定义与此处一致，之后不再赘述。对于同一种算例来说，其升频扫描得到的幅频响应特性曲线与降频扫描得到的幅频响应特性曲线在两边的部分(即远离共振区部分)应该是重合的。根据这个特点可从图中区分各条曲线。

由图5可知，当激波振荡频率接近主共振频率时，系统振动幅值急剧增大，同时，在共振峰处，正向扫频和反向扫频过程均存在明显的振动突跳和双稳态现象，在升频过程突跳点A与降频过程突跳点B之间，形成双稳态区间。从图5(a)，随着振荡激波的强度不断增大，系统共振峰幅值不断增大，主共振幅频响应曲线不断整体右移，主共振频率不断增大。对于弱激波而言，σ=17°，系统不存在振动突跳和双稳态现象；当激波强度增大即σ=18°时，系统存在明显的振动突跳和双稳态现象，共振峰处的双稳态区间对应的频率为ω∈[10.3,13.0]；当σ=20°时，升频跳跃点较之前右移，而降频跳跃点大幅右移，使得双稳态区间变窄，此时共振峰处的双稳态区间对应的频率为ω∈[14.2,15.1]；继续增大激波强度，当σ=22°，主共振幅频响应曲线较之前整体右移，但振动跳跃和双稳态现象再次消失；当进一步增大激波强度，考虑强激波σ=24°时，系统出现显著的振动突跳和双稳态现象，不仅共振峰值和主共振频率显著提高，双稳态区间宽度大幅变宽，此时共振峰处的双稳态区间对应的频率为ω∈[19.9,25.1]。

(a)

3.2 激波振荡幅值对受热壁板主共振特性的影响

分别取激波振荡幅值为A=0.005、0.010、0.030、0.050、0.100、0.200，来流马赫数为Ma=3.5，激波角为σ=18°，温升为ΔT/ΔTcr=0.10，中心位置为ξ0=0.50，绘制系统的幅频响应特性曲线如图6所示。

由图6可知，随着激波振荡幅值的增大，系统共振峰值不断增大，并且发生上升跳跃和下降跳跃的频率均随着激波振荡幅值的增大而增大。当激波振荡幅值较小时，即A=0.005, 系统不存在振动突跳和双稳态现象；当激波振荡幅值增大时，即A=0.01时，系统出现明显的振动突跳和双稳态现象，共振峰处的双稳态区间对应的频率为ω∈[11.1,13.0]；当A=0.03时，升频跳跃点和降频跳跃点都像右移，双稳态区间变宽，此时共振峰处的双稳态区间对应的频率为ω∈[13.4,16.1]；当激波振荡幅值进一步增大时，即A=0.050、0.10，从图6还可知，随着激波振荡幅值的增大，升跳跃点和降频跳跃点不断右移，但双稳态区间不断减小直至消失。当A=0.050时，共振峰处的双稳态区间对应的频率为ω∈[16.0,16.9]；当A=0.100时，振动突跳和双稳态现象消失。而当继续增大激波振动幅值，即A=0.200时，从图中可以看出，系统出现显著的振动突跳和双稳态现象，并且双稳态区间宽度将急剧变宽，此时共振峰处的双稳态区间对应的频率为ω∈[20.1,29.4]。

图6 不同激波振荡幅值下受热壁板的幅频响应曲线

3.3 温度对振荡激波作用下受热壁板主共振特性的影响

分别取无量纲温升为ΔT/ΔTcr=0、0.1、0.5、1.0，来流马赫数为Ma=3.5，激波角为σ=18°，中心位置为ξ0=0.500，激波振荡幅值为A=0.010，绘制系统的幅频响应特性曲线如图7所示。

从图7可知，随着温升的不断增大，系统共振峰值不断增大。当不考虑温升时，即ΔT/ΔTcr=0，系统不存在振动突跳和双稳态现象；当考虑温升时，即ΔT/ΔTcr=0.1、0.5、1.0时，系统存在明显的振动突跳和双稳态现象；从图7还可知，双稳态区间的宽度随着温升的增大几乎不改变，发生上升跳跃和下降跳跃的频率随着温度的不断增大而减小，系统双稳态区间位置不断向左移，即系统主共振区对应的频率逐渐减小，这说明了提高温度对系统产生“刚度弱化效应”。

(a)

3.4 激波振荡的中心位置对受热壁板主共振特性的影响

分别取激波振荡的无量纲中心位置为ξ0=0.100、0.300、0.45、0.500、0.700、0.900，来流马赫数为Ma=3.5，激波角为σ=18°，温升为ΔT/ΔTcr=0.1，激波振荡幅值为A=0.010，绘制系统的幅频响应特性曲线如图8所示。

图8 不同中心位置下受热壁板的幅频响应曲线

从图8可知，激波振荡的中心位置对系统共振峰值有很大的影响，当激波围绕壁板中点振荡时，即ξ0=0.50时，系统共振峰值最大。随着中心位置向两端移动时，系统共振峰值不断减小。相反地，随着中心位置向两端移动时，主共振频率不断增大。当激波围绕壁板中点振荡时，即ξ0=0.50时，系统出现明显的振动突跳和双稳态现象；而当中心位置移动到ξ0=0.45时，双稳态区间的宽度急剧变窄；当中心位置移动到ξ0=0.1、0.9、0.3、0.7时，振动突跳和双稳态现象消失。由此可见，将激波振荡中心位置往两端移动时可有效地抑制振动突跳和双稳态现象。

3.5 来流马赫数对受热壁板主共振特性的影响

分别取来流马赫数为Ma=3.0、3.5、3.8、4.0、4.5，来流动压为λu,l=30，激波角为σ=20°，温升为ΔT/ΔTcr=0.1，中心位置为ξ0=0.500，激波振荡幅值为A=0.01，绘制系统的幅频响应特性曲线如图9所示。

由图9可知，随着来流马赫数不断增大，系统共振峰幅值不断增大，主共振幅频响应曲线不断整体右移，主共振频率不断增大。从图9可知，来流马赫数对受热壁板主共振特性的影响与振荡激波强度、激波振荡幅值对受热壁板主共振特性的影响类似。当来流马赫数较小时，Ma=3.0，系统不存在振动突跳和双稳态现象；当马赫数增大，即Ma=3.5，系统存在明显的振动突跳和双稳态现象；当来流马赫数增大，升频跳跃点较之前右移，而降频跳跃点大幅右移，双稳态区间宽度不断变窄直至消失；而来流马赫数较大时，即Ma=4.5，系统出现显著的振动突跳和双稳态现象，不仅共振峰值和主共振频率显著提高，双稳态区间宽度大幅变宽。

图9 不同来流马赫数下受热壁板的幅频响应曲线

4 结 论

本文基于Von-Karman 大变形理论和当地一阶活塞流理论, 建立了振荡激波作用下受热壁板振动的动力学模型，通过4阶龙格-库塔法对非线性动力学模型进行求解，得到了振荡激波作用下受热壁板的幅频特性响应曲线，并详细讨论了振荡激波强度、激波振荡幅值、温度、激波振荡的中心位置、来流马赫数对系统振动突跳和双稳态特性的影响。主要结论如下：

(1)受热壁板在振荡激波作用下存在振动突跳和双稳态现象等典型的非线性动力学行为，并且系统的幅频响应特性曲线在双稳态区表现出明显的“硬特性”。

(2)随着振荡激波强度、激波振荡幅值和来流马赫数增大，系统共振峰值不断单调增大，而双稳态区间由不存在到逐渐变大，再变小直至消失，然后再急剧变大。

(3)温度增大会对系统产生“刚度弱化效应”，使主共振幅频响应曲线整体左移，但对系统的双稳态区间影响较小。随着温度的增大，系统共振峰值不断增大。

(4)当激波围绕壁板中点振荡时，系统共振峰值最大。随着中心位置向两端移动时，系统共振峰值不断减小，而主共振频率不断增大。将激波振荡中心位置往两端移动时可有效地抑制振动突跳和双稳态现象。本文主要是从数值的角度研究了振荡激波作用下受热壁板的主共振特性，在我们下一步的研究工作中，将从试验的角度研究受热壁板在振荡激波流场中的主共振特性。