安效民，冯家悦，周 悦，孙 伟

(西北工业大学航天学院空天飞行技术研究所，西安 710072)

超音速飞行器的翼、身及重复使用火箭的箭体多采用轻质薄壁加筋结构，其一侧直接承受气动载荷作用，另一侧为腔体结构。这种薄壁结构在气动载荷作用下发生变形，导致外部流场的边界发生改变，引起流场的结构(如边界层、激波等)和气流参数(如速度、压强等)的改变，使得作用在结构上的气动载荷发生变化，进而使结构产生新的变形或者振动。这种结构弹性、气动载荷与惯性载荷间的相互耦合，会导致壁板气动弹性问题，出现不稳定性现象[1-4]。

与传统翼、舵等升力面构型不同的是，薄壁结构在跨音速流表现出了双重的非线性流固耦合特征。一方面，壁板位移响应通常与壁板厚度为相同量级，壁板面内应力存在，弯曲和拉伸之间会发生耦合，在周围受约束(固支或简支)条件下，其形变一般不会引发结构的迅速破坏，但会引起结构疲劳，表现为几何非线性[5-7]。另一方面，壁板的响应还受跨音速气动非线性的影响：① 由于壁板向上或向下的振荡，会增强或减弱激波强度，并且使得激波前后运动，激波运动可能是连续的，也可能是间歇的，或者连续和间歇互相转换，从而导致壁板复杂的响应形态在相对较宽的动压范围内持续存在，而且会有多种形态之间的演化；② 激波与边界层之间的干扰可能造成流动分离，激波或分离涡的运动进一步加剧了流动的动态非线性特征，使得壁板响应更为复杂[8-10]，而且在跨音速区，粘性效应本身对壁板响应具有增稳或失稳作用[11-13]。

在这些跨音速流中强非线性因素的作用下：一方面壁板的稳定边界呈现出与流线型升力面构型(如翼、舵部件等)相似的跨音速凹坑[5]；另一方面，在进入到不稳定区域后，壁板表现为多个不稳定屈曲、极限环颤振或者更为复杂的振荡行为，诸如周期性、拟周期、非周期和混沌等复杂响应[14]。准确预测和确定飞行器在跨音速下壁板结构的稳定边界，并分析其在不稳定域内的形态演化规律，有助于揭示壁板各类复杂行为的诱发机理，寻求抑制甚至消除最严重形态响应的方法，可有效降低飞行器的疲劳损伤，尤其是对于可重复使用火箭，其在上升段有较长时间飞行在跨音速及低超音速范围，薄壁结构受到气动载荷作用而产生的变形或振动对结构疲劳有着重要影响。

对于跨音速流固耦合效应下壁板的非线性响应问题，当前无法建立统一的动力学模型，传统基于变分方法或其他离散化方法，将偏微分方程化为常微分方程组，进而进行特征值性态分析的应用受限。DOWELL[15]提出了4 类壁板颤振分析理论，后来CHENG 和MEI[16]将其扩充为6 类，他们都指出：在跨音速范围内，使用非线性结构理论和求解流体的Euler/Navier-Stokes 方程是一种有效途径。

DAVIS 等[8]联立求解Euler 方程和非线性板的结构方程，研究了M∞=0.8～2.5 范围内二维壁板的气弹响应，发现激波的出现会导致气弹响应出现发散和极限环等现象。GORDNIER 等[10]的研究中利用可压N-S 方程求解壁板表面的气动力，并结合Von-Karman 理论，采用隐式迭代求解壁板响应，分析了M∞=0.8下壁板的稳定性，并解释了粘性边界层效应对气动弹性稳定性和颤振失稳后形态的影响。HASHIMOTO 等[12]的研究则利用了耦合求解思路，将N-S 方程和Von-Karman 方程联立起来，研究了边界层效应对壁板颤振的影响。ALDER[13]的研究中采用隐式有限体积法对N-S 方程进行了求解，耦合考虑了几何非线性的壁板有限元模型，分析了从高亚音速到低超音速下湍流边界层对壁板系统稳定边界的影响。SHISHAEVA等[17-18]联立ABAQUS 结构求解器和Flowvision流场求解器，分析了二维壁板从高亚音速到低超音速阶段的壁板响应，指出随马赫数变化中，壁板会呈现为屈曲、动态稳定、单模态颤振、耦合颤振等多种响应形式，并分析了加速、减速效应下壁板的响应形态变化。

国内肖艳平等[19]基于一阶活塞气动力理论，采用伽辽金法分析了边界松弛对超音速气流壁板颤振响应的影响，结果表明：随边界约束的松弛，壁板可能产生颤振极限环振动。吴志强等[20]通过数值计算庞加莱映射分岔的方法，讨论了翼型在不可压流中的极限环颤振随气流速度变化引起的分岔行为，给出了8 种典型的相图和谱图，并分析了闭轨分岔的诱因。钮耀斌等[21]采用椭圆函数谐波平衡法研究分析了翼型的超音速非线性颤振问题，结果表明极限环振荡临界速度随着弹性轴位置与翼弦中点距离的减小而不断增大，且随着重心位置与弹性轴距离的增大，极限环振荡临界速度存在一个极小值点。可以看出，这些非线性颤振大多关心超音速或低速流动特点，气动力建模考虑了线性模型。朱世权等[22]基于CFD/CSD 单向流固耦合计算方法，对0.4～1.6 之间8 种马赫数下大展弦比机翼进行了静气动弹性数值研究，结果表明：机翼翼尖位移和机翼最大应力在跨音速范围发生了突变。可为相关大展弦比机翼的设计与分析提供参考。刘燚等[23]采用曲面涡格法对柔性飞机进行非线性静气动弹性分析，结果表明：曲面涡格法在可压缩情况下载荷计算精度较好且气动力曲面建模优势明显，可用于工程复杂模型的曲面气动力计算。

本文针对壁板在跨音速流中的稳定性和流固耦合形态演变规律展开研究，分析壁板在不同结构和流动参数下随马赫数变化过程中的形态演化规律；考虑非定常加速效应的影响，分析不同加速比下的形态演化；考虑粘性效应的影响，分析不同边界层厚度下的形态演化。

1 基于CFD/CSD 耦合的壁板非线性响应分析方法

为了求解由气动和结构两场非线性引起的非线性壁板非线性响应问题，采用了CFD/CSD 时域耦合求解方法[24-25]，其基本流程(图1)如下：

图1 CFD/CSD 耦合流程Fig. 1 Sketch of CFD/CSD coupling method

式中：S为流场和结构系统之间的转换矩阵[26]；x(n+1/2)为在第 (n+1/2)时刻流体域中表面网格节点的位移。

2) 利用动网格技术更新流体网格。考虑惯性力影响，在流固耦合边界上设置压力梯度条件：

2 壁板模型与算法验证

2.1 壁板模型

二维壁板几何模型和计算域如图2 所示，壁板长度为a、厚度为h。弹性薄壁板边缘固定并与刚性表面平滑连接，薄壁板的边界条件为固支。当在无粘条件下计算时，刚性表面设置为无穿透条件，L1=L2=L3=10a。当考虑粘性条件时，刚性表面设置为无滑移条件，并且使L1具有恰当的长度，使求解薄壁板处的流场时，可以形成具有所需厚度的粘性边界层，如图3 所示。

图2 壁板的几何模型和计算域Fig. 2 Panel model and computational domain

图3 壁板粘性边界层的形成Fig. 3 Viscous boundary layer of the panel

来流马赫数和密度定义为M∞和 ρ∞。壁板基准模型参数设置为[17]：a=0.3m，a/h=300，E=2×1011Pa ，泊松比ν=0.3 ，壁板密度ρs=7800 kg/m3。为方便对比分析，定义如下参数：

2.2 壁板气动弹性响应验证

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

研究了不同网格尺寸和时间步长下M∞=1.12，μ=1.64×10-4时的壁板响应。壁板结构的有限元模型由20 个单元组成，考虑了壁板的两端固支边界条件。计算工况如表1 所示。

表1 网格数量和时间步长的计算工况Table 1 Calculation conditions of grid size and time-step

图4 显示了在不同网格数量下的计算结果在壁板参考点(x/a=0.75)处的时间历程曲线，可以看到，在不同网格数量下，壁板极限环振幅完全相同的，随着时间推移，瞬时相位有微小的差异。图5 显示了在不同时间步长下的计算结果，发现在不同时间步长下，极限环幅值完全相同，但随着时间推移，极限环也有微小相移。综合考虑计算精度和效率，本文最后选用的网格数量为161×101，时间步长为0.00001 s。

图4 不同网格数量下的壁板响应Fig. 4 Response of panel under different number of grids

图5 不同时间步长下获得的壁板响应Fig. 5 Response of panel under different time-steps

2.2.2 二维壁板跨音速稳定性分析验证

为验证本文算法的可靠性，壁板参数与文献[6,8,10,13]相同：a/h=50，E=7×1010Pa，泊松比ν=0.3 ，质量比μs=0.1。

计算了马赫数从0.8～1.3 时壁板的稳定边界，并引用了文献[6,8,13]的结果，绘制在图6 中进行比较。如图所示，当M∞≤1时，壁板的不稳定形式表现为屈曲，当M∞＞1时，壁板不稳定表现为颤振及极限环振荡等形态。这种不稳定的边界结合到一起，也会表现出一般升力面构型(如翼、舵结构等)在跨音速流中的颤振凹坑现象。图6 显示了土星五号的飞行轨迹，可以看出，在M∞≈1的区域内，不稳定的临界动压明显低于飞行动压。DOWELL[6]、DAVIS 等[8]及ALDER[13]的研究都表明了这一点，从计算结果来看，本文与DAVIS 和ALDER 等的结果相符较好，由于DOWELL 的计算中采用了势流理论(在M∞=1时有奇异)，有较大差异。图7显示了M∞=0.9时不同动压下参考点的屈曲位置，可以看出，在不同的扰动下，壁板屈曲的平衡位置有两个。当动压较小时，这两个屈曲位置相对于壁板初始位置是对称的，随着动压增大，屈曲位置呈现出不对称，这种不对称与较大动压下不同扰动所造成的激波强弱有关，从该算例在跨音速范围内的计算来看，本文所使用的方法与文献[10]对比符合较好。

图6 壁板在跨音速范围内的稳定边界Fig. 6 Stable boundary of panel in transonic domain

图7 壁板参考点位置随动压的变化Fig. 7 Deflection of the panel vs. dynamic pressure

3 跨音速流中定常状态下壁板流固耦合形态演化分析

3.1 计算条件

考虑了无粘条件下的3 种壁板模型，三种计算工况如表2 所示。基准状态壁板模型参数设置如第2.1 节所示，而针对壁板长厚比和来流密度(来流动压)变化的壁板模型见表2。初始扰动为(考虑了结构前2 阶模态形式)：v=0.001sin(πx/a)+0.001sin(2πx/a) 。 马赫数范围为 0.7≤M∞≤2.0。利用CFD/CSD 耦合方法计算壁板响应，通过壁板参考点(x/a=0.75)的响应历程、频谱图、庞加莱映射图、相平面图等分析和判别壁板形态。

表2 三种工况的参数Table 2 Calculation conditions of three panel models

3.2 定常情况下壁板的形态演化

壁板的形态随着马赫数增大的变化过程如图9 所示。图10 和图11 分别显示了参考点变形位置/振荡幅值、振荡频谱峰主导频率随马赫数的变化，可以发现：

图9 随马赫数变化下三种壁板的形态演化Fig. 9 Morphological evolution of three kinds of panels with Mach number

图10 变形/幅值随马赫数的变化Fig. 10 Deflection/ Amplitude of the panel vs. Mach number

图11 振荡主频随马赫数变化Fig. 11 Main frequency of oscillation vs. Mach number

1) 对于三种壁板计算工况，随马赫数的逐渐增大，壁板的响应形态大体上从稳态收敛、第一模态极限环(LCO)振荡、屈曲、跨音速颤振、非共振型极限环(LCO)、共振型极限环 (LCO)、高频周期振荡、高频非周期振荡到稳态收敛的过程，这与文献[17]的结果大体一致。

2) 第一模态极限环振荡、跨音速颤振和非共振型极限环振荡、共振型LCO 的形态如图12～图14 所示，可以观察到，第一模态极限环表面为纯单模态振荡形式，其振荡频率为13 Hz，远远小于一阶固有频率ω0=61.76 Hz。跨音速颤振会表现出明显的非对称延迟现象，这是由流场中存在的多激波结构，以及激波前后移动所导致(如图15所示)，此时频谱由多个峰值组成，主频在壁板固有频率附近，这种现象在GORDNIER 等[10]和ALDER[13]的研究中也有所体现。当壁板厚度增加(Case 2)、或者来流密度变小(Case 3)，并未观察到明显带有延迟振荡的跨音速颤振。随马赫数增加，延迟振荡变弱，呈现为如图14 所示对称形式，振荡频谱的3 个峰值近似为1∶2∶3 组成，其中第二模态是一个小的分量，随着马赫数的增加，该分量比例逐渐增加，向共振型LCO 过渡。其演化机理大体如下：随着马赫数的增大，气流的能量也增大，气流能量向结构传递的过程中，可能会引起其第一模态的增长，第一模态增长过程中，振荡频率也会增加，存在一个临界值，之后会激发第二模态，其增长通过与第一模态的共振来维持，随着气流能量的继续增加，第二模态的主导作用也越来越明显，当流入壁板的能量等于流出壁板的能量时，第一、二阶模态的振幅达到稳定，进入极限环振荡。相似的分析可见文献[27]。

图12 工况1 第一模态LCO 的响应和频谱图Fig. 12 Time-histories and frequency spectrum of first-mode LCO (Case 1)

图13 工况1 跨音速颤振的响应和相图Fig. 13 Time-histories and phase portraits of transonic flutter (Case 1)

图14 工况1 非共振型LCOFig. 14 Time-histories and phase portraits of nonresonant LCO (Case 1)

图15 工况1 跨音速颤振瞬时φ=180°的流场压强Fig. 15 Pressure distributions at φ=180° of transonic flutter (Case 1)

3) 亚音速下的第一模态LCO 和跨音速颤振中，其振荡主频低于壁板第一阶固有频率值，表现为单模颤振，与经典颤振中两个或者多个结构模态耦合形成机理不一样的是，其形成是由于气流的不稳定引起的，因此产生的颤振模态可能同时来源于壁板的结构模态和气流模态。当壁板发生共振极限环发生时，表现为耦合模态振荡，如图16 所示，此时随马赫数的增加，主频增长减慢(图11)。

图16 工况1 共振型LCO 的响应和频谱图Fig. 16 Time-histories and frequency spectrum of resonant LCO (Case 1)

4) 对于高频周期和非周期振荡而言，振荡中存在多个高频模态，此区域很可能是疲劳损伤最严重的区域：即使偏转幅值与第一模态极限环相近，但其壁板形态包含较高的模态，因此应力振幅比第一模态极限环的应力振幅高得多，此外，该区域包含较高的频率，因此一旦产生，极有可能会在短时间内迅速积累大量的疲劳损伤。

5) 在更大的壁板厚度(工况2)和更小的动压(工况3)下，壁板演化过程与基准状态(工况1)相比主要有3 点不同：

a) 在亚音速区没有第一模态极限环振荡；

b) 由于壁板较厚(或由于无量纲气流密度较小，即动压较小)，M∞=1时无法激发出跨音速颤振，反而保持稳定状态；

c) 在低超音速区，作为过渡的高频周期振荡和第一模态极限环振荡没有激发，表现为高频非周期振荡。

3.3 考虑非定常加速效应下壁板的形态演化分析

3.3.1 计算条件

针对三个计算工况，定义非定常加速中马赫数的变化如下[18]：

其中，M1=0.7 和M2=1.7。通过调整时间T来改变流动的加速度。为了使得壁板响应进入M1和M2之间前得到较为稳定的形态，选取适当的M0即可。分析中取T=10 s、5 s 和2.5 s，代表了加速度从小到大。需要注意的是，该加速条件是通过远场边界实施，加速效应对壁板产生的影响会存在一定的时间延迟，但总体延迟时间很小，分析的几种工况不大于0.004 s。

3.3.2 加速条件下壁板的形态演化

当考虑非定常加速效应时，壁板的形态演化如图17 所示；参考点偏转的时间历程如图18 所示；庞加莱映射如图19 所示，从图中可以清楚地看到壁板块动力学的大部分分岔。不同加速情况下振幅和主导频率与马赫数的关系比较如图20 所示(灰色点线图为定常马赫数下壁板响应的幅值和频率)，可以看出：

图17 不同加速度下壁板的形态演化(工况1)Fig. 17 Morphological evolution of panels under different accelerations (Case 1)

图18 不同加速情况下的壁板响应历程Fig. 18 Response history of panel under different acceleration

图19 不同加速情况下的庞加莱映射Fig. 19 Poincare maps under different accelerations

1) 在亚音速情况下，所有加速度下的响应幅值都接近于定常流中的幅值，加速度越大，屈曲区的值越大，且马赫数越接近1，不同加速度下幅值差异越小(图20(a))。

图20 不同加速度条件下的振幅和主导频率Fig. 20 Amplitude and dominant frequency under different accelerations

2) 在跨音速时，具有非定常加速效应时发生振荡所需的马赫数略大，且加速度越大，振荡发生的马赫数越大，即振荡延迟越明显。

3) 在 1≤M∞≤1.26时，具有非定常加速效应和定常时的结果吻合的很好，不仅振幅接近，而且主导频率只有第一和第二频率。

4) 对于更高的马赫数，存在明显差异：

a) 不存在作为过渡的高频周期振荡；

b) 由于更高频率的振荡需要较多的时间去发展，当加速度较低时(T=10 s、T=5 s)，高频振荡发生的马赫数大于定常流发生高频振荡所对应的马赫数，振荡幅值小于定常流，且随着加速度的提高，高频振荡发生的马赫数逐渐提高，振荡幅值减小，但高频区频谱中高频占比要大于定常流流，且频率值随马赫数增加有一个明显的先快速提高后快速下降的过程(图20(b))；

c) 对于较高的加速度(T=2.5 s)，较高频率发展所需的时间显得有些过长，在高频分量明显显现之前，流速已经跨越可以支撑高频分量继续发展的马赫数域，因此高频振荡将不会产生；

d) 当T=5 s 时，存在一个过渡状态—高频准周期运动，其频率主导模态为六阶模态；

e) 当T=2.5 s 时，共振极限环会演化为一个第一、二阶模态交替的新形态，演化较为平滑，无法准确确定分岔点，在该阶段壁板形状有两个形态，分别接近第一阶振型和第二阶振型，两种形态交替出现交替主导。当M∞≈1.7 时，振荡幅值极小，接近消失。

图21 和图22 显示了壁板厚度增加(工况2)和来流密度减小后(工况3)的非定常计算结果，可以看出：

图21 不同加速度下壁板的形态演化(工况2)Fig. 21 Morphological evolution of panels under different accelerations (Case 2)

图22 不同加速度下壁板的形态演化(工况3)Fig. 22 Morphological evolution of panels under different accelerations (Case 3)

1) 在亚音速、跨音速和马赫数较低的低超音速域，可以得到与基准状态相同的结论。

2) 在马赫数较大时，对于工况2，在加速情况下，当M∞≥1.38时壁板幅值将缓慢减小，直到马赫数接近1.7 才会恢复为屈曲状态。当加速度较小(T=10 s)时，在这一区域将存在一个频率略小于定常流、幅值初始值接近最大幅值且随马赫数增加而逐渐减小的高频非周期运动；当加速度较大(T=2.5 s)时，高频区被跨过，相应的马赫数域由一个幅值随马赫数增加而逐渐减小的第一模态极限环振荡所代替。

3) 在马赫数较大时，对于工况3，当加速度较小(T=10 s)时，存在一个马赫数域较广(类似于壁板更厚的情况)但幅值较低(类似于基准状态)、频率较低的高频非周期振荡；当加速度较大(T=2.5 s)时，高频区被跨过。

3.4 考虑粘性效应下壁板形态的演化分析

3.4.1 计算条件

为了分析粘性效应对流固耦合形态演化规律的影响，取壁板中心点距上方99%远场速度处位置为附面层的外边界。本文针对两种不同附面层厚度 δ/a=0.025 和 δ/a=0.05进行分析，为取得对应的附面层厚度，通过调整壁板前端到远场的距离L1，来发展形成附面层，如图3 所示。

3.4.2 定常条件下壁板的形态演化

定常条件不同附面层厚度下壁板的形态演化如图23 所示，振幅和主导频率与马赫数的关系比较如图24 所示(灰色点线图为不考虑粘性时壁板的幅值和频率)，发现当考虑粘性后：

图23 不同附面层厚度下壁板的形态演化Fig. 23 Morphological evolution of panels under different boundary layer thickness

图24 不同附面层厚度下的振幅和主导频率Fig. 24 Amplitude and dominant frequency under different boundary layer thickness

1) 亚音速下的稳态和第一模态LCO 转变为屈曲，跨音速范围内部分马赫数下的动压不足以支撑振荡的产生，壁板的振荡恢复为屈曲状态。

2) 较高的模态会产生较大的阻尼，使得在无粘模型高频区用来使较高模态发展增长的这一部分能量被耗散掉，高频区范围大大缩小，当附面层为 δ/a=0.025时，高频区完全消失，壁板形态由共振极限环转换为低频准周期振荡(主导模态仍为第一和第二模态，但存在多个峰值)和第一模态极限环(壁板运动接近一阶振型)；随着附面层的继续增厚 δ/a=0.05时，高阶模态的影响将完全消失。

3) 随着附面层的增厚，壁板振荡的最大幅值和频率都会有所减小，而在较高马赫数时(M∞≥1.5)，幅值和频率反而有所增大。

3.4.3 加速条件下壁板的形态演化

考虑非定常加速效应后(T=10 s、T=2.5 s)不同附面层厚度下壁板的形态演化如图25 和图26所示，振幅和主导频率与马赫数的关系比较如图27和图28 所示，结果表明：

图25 不同加速度下壁板的形态演化(δ/a=0.025)Fig. 25 Morphological evolution of panels under different accelerations (δ/a=0.025)

图26 不同加速度下壁板的形态演化(δ/a=0.05)Fig. 26 Morphological evolution of panels under different accelerations (δ/a=0.05)

图27 δ/a=0.025时不同加减速条件下的振幅和主导频率Fig. 27 Amplitude and dominant frequency under different accelerations ( δ/a=0.025)

1) 考虑非定常加速效应后，振荡起始和结束的马赫数会向较大的马赫数偏移，附面层越厚偏移越明显。

2) 随着加速度的提高，最大幅值和最大频率均有所减小，附面层越厚越明显。

3) 当 δ/a=0.025时，在原高频域，当具有加速效应后，由于附面层厚度不足导致的演化不再出现。

4) 考虑非定常效应后振荡的主频随马赫数呈现出台阶式变化，即在一定马赫数范围内，振荡主频保持不变，加速度越高、附面层越厚，这种台阶状变化越明显(图27(b)、图28(b))，这与考虑粘性效应后阻尼增大有关。

图28 δ/a=0.05时不同加速条件下的振幅和主导频率Fig. 28 Amplitude and dominant frequency under different accelerations ( δ/a=0.05)

4 结论

本文基于CFD/CSD 耦合方法，分析了跨音速范围内壁板形态的演化规律。从仿真结果分析总结出如下规律：

(1) 在跨音速范围内，随着马赫数的增加，壁板的响应呈现出复杂的演化形态：稳态收敛、第一模态LCO、屈曲、跨音速颤振、非共振型LCO、共振型LCO、高频周期振荡、高频非周期振荡等，其中亚音速下的第一模态LCO、跨音速颤振和非共振型LCO 的振荡主频小于结构1 阶固有频率，呈现出单模态颤振的形态。

(2) 从壁板随马赫数的变化来看，非周期振荡区域出现了高振幅和高频率特征，很可能是典型的单模颤振中疲劳损伤最严重的区域。

(3) 考虑非定常加速效应后，各演化阶段向马赫数更大的方向推移。当来流加速度增加时，发生高频振荡的马赫数范围会变短，如果加速度足够大，则高频振荡形态会消失。

(4) 当考虑粘性效应后，高频区范围将大大缩小，并随着附面层的增厚，高频区最终完全消失。且振荡的主频呈现出台阶式变化，加速度越高、附面层越厚，这种台阶状变化越明显。

(5) 总体来看，当考虑非定常加速效应和粘性效应后，对于单模态颤振发生会延迟，对于引起壁板疲劳损伤最为严重的高频振荡会有所抑制，这在飞行器设计中是有利的因素，传统基于无粘假设和定常分析的结果偏于保守。