基于流固耦合的激波/边界层干扰作用下壁板颤振特性

2023-04-19刘为佳李映坤陈雄李春雷

航空学报 2023年6期

刘为佳，李映坤，陈雄，李春雷

南京理工大学机械工程学院，南京 210094

壁板颤振和激波/边界层干扰是超声速飞行器结构设计与优化的重要基础科学问题。壁板颤振是飞行器表面薄壁结构在惯性力、弹性力和气动力耦合作用下发生的一种自激振动[1］，一般呈现出有限幅值的极限环颤振，其非线性动力学特性会造成结构的疲劳损伤累积，影响壁板结构的疲劳寿命，甚至对飞行器的飞行性能和安全产生不利影响。作为经典的气动弹性问题，Dowell[2-3］、Mei[4］和杨智春[5］等学者全面总结了经典壁板颤振研究所采用的气动力分析模型、结构分析模型、降阶模型、数值求解方法及理论分析和实验方面的研究成果。

激波/边界层干扰是超声速/高超声速飞行器内、外流场中普遍存在的一种复杂流动现象，既会急剧增强壁面的湍流脉动和黏性耗散，引发升力降低、阻力增加、失速提前等不良后果，又会导致进气道喉道附近流场畸变、总压恢复下降、发动机不起动等[6］。对于激波/边界层干扰现象，众多学者展开了大量研究，涉及流场结构、物理机制、流动控制方法、低频不稳定性、侧壁三维效应等[7-13］。同时，目前超声速/高超声速飞行器一般采用薄壁结构以减轻重量，在发动机流道内部，流场中存在激波、膨胀波、激波反射与相交等复杂波系，而在飞行器外流环境下，不同部件之间会存在多体干扰，一个部件的激波会作用在其他部件表面。因此，对于超声速/高超声速飞行器，不可避免的存在激波/边界层干扰流动与弹性壁板结构之间的流固耦合现象。

在激波/边界层干扰主导流动与弹性壁板相互作用的实验方面，学者做了大量研究。Daub等[14］采用电容式位移传感器测量壁板变形振动，指出激波冲击位置快速变化下壁板的最大振幅约1 mm；Varigonda 等[15］发现当壁板底部为环境压力时弹性壁板呈现凸起变形，而当壁板底部为来流静压时弹性壁板呈凹曲变形并会引起局部流动分离；Tripathi 等[16-17］重点研究了来流雷诺数和激波冲击位置对壁板动力学响应的影响；Spottswood 等[18］则提出了基于快速响应压敏涂料和激光测振的全域非接触式测量技术，并获取了激波/边界层干扰流动与弹性壁板耦合系统的压力脉动和结构振动响应；Tan 等[19］采用多个激振器驱动弹性壁板连续变形，分析了壁板曲率和激波冲击位置对压力分布和边界层分离区长度的影响，利用弹性壁板的弯曲变形可对激波/边界层干扰流动进行控制。

综上所述，本文采用有限体积法求解非定常可压缩Navier-Stokes 方程，获得了壁面的气动力载荷，用有限差分法求解基于几何大变形理论的弹性壁板运动方程，并采用双向流固耦合计算方法对激波/边界层干扰下壁板的气动弹性响应过程进行数值模拟，分析了不同激波冲击位置下弹性壁板的颤振特性以及激波/边界层干扰流动分离特性，研究为激波主导流动中弹性壁板结构设计及激波/边界层干扰流动分离的控制提供理论基础。

1 控制方程与计算方法

1.1 流体区域

针对激波/边界层干扰作用下弹性壁板的气动弹性数值模拟研究，Visbal[20-21］采用双向耦合方法，研究了斜激波主导流动下弹性壁板的动力学特性，发现当入射激波强度足够强时，弹性壁板发生颤振所需的临界动压更低，并分析了壁板振动对流场分离区的影响；Boyer 等[22-23］将Visbal 的计算模型拓展到三维，分析了增压比和动压对壁板颤振及激波/边界层流动分离的影响；Shinde等[24-25］采用基于大涡模拟方法研究了变形壁面对激波/边界层干扰的控制作用；Bhatia 等[26］研究了壁板曲率引起的气动载荷非线性特征对跨声速流动中壁板的颤振影响；Whalen 等[27］研究了高速压缩拐角流动中激波/边界层干扰下弹性斜坡的流固耦合作用；Li 等[28］研究了壁板反馈控制速度和边界层厚度对壁板颤振和流动分离特性的影响；An[29］和李映坤[30］等采用流固耦合方法，研究了二维曲壁板在斜激波冲击作用下的非线性气动弹性。 Wang 等[31］采用HODMO（Higher-Order Dynamics Modes Decomposition）分解方法对跨声速流动下壁板的模态稳定性进行了分析。然而，针对激波/边界层干扰下壁板的气动弹性问题，激波冲击位置对壁板振动响应特性以及边界层流动分离特性的影响尚未见公开报道。

激波/边界层干扰作用下弹性壁板振动会引起流体计算区域的不断变化，流动控制方程需要考虑计算区域网格的运动。本文采用基于ALE（Arbitrary Lagrangian Eulerian）方法的可压缩非定常Navier-Stokes 方程来描述，即

式中：Q 为守恒变量；Fc为无黏通量；Gv为黏性通量；n 为控制体表面法向量，以上各式的具体形式参考文献[32］。

基于多块结构网限体积法进行流动控制方程的离散，对流通量计算采用三阶MUSCL（Monotone Upstream centered Schemes for Conservation Laws）重构和AUSMPW+（Advection Upstream Splitting Method by Pressurebased Weight functions）格式[33］，黏性项计算采用二阶中心差分格式，时间推进采用双时间步LU-SGS（Lower Upper Symmetric Guass Seidel）离散方法[34］。

1.2 固体结构区域

由于壁板振动位移较大，本文采用考虑几何非线性二维弹性壁板运动方程[35］描述，即

式中：X(x，y)为壁板上某一位置的坐标；h、ρs和l表示壁板的厚度、密度和弧长；n为单位法向量；t=∂X∂l为单位切向量；p为作用在壁板表面的气动力，其中表现几何非线性项的横向应力q为

式中：M为弯矩；κ为当地曲率；EB为弯曲模量。

壁板两端设置为简支边界条件，即

描述壁板几何非线性项的面内张力τ为

式中：Es为壁板的拉伸系数；l0为壁板初始弧长。

采用二阶有限差分格式对式（2）中的空间导数离散，时间推进采用欧拉法。

1.3 双向流固耦合方法

壁板颤振过程中的双向流固耦合问题采用分区迭代耦合算法求解，如图1 所示。该算法是指在每一时间步对流体和固体单物理场进行一次求解，并通过将流体域气动压强传递到固体域，而固体域将振动变形的位移和速度量传递到流体域，在时间步推进过程中实现耦合问题的求解。本文计算模型的流固耦合界面网格相匹配，无需进行特殊处理，另外壁板颤振过程中运动位移较小，流体区域网格运动只需根据运动后的边界进行调整[21］，表达式为

图1 松耦合算法求解时序图Fig.1 Generic cycle of loosely coupled algorithm

式中：δwi，j为壁板表面的振动位移；s为沿计算坐标的弧长，网格在jmax之后保持不变形。

1.4 数值计算方法验证

为验证数值方法的准确性，本文研究了来流压比p3/p1=1.4、无量纲动压λ=875、来流马赫数Ma=2.0 条件下斜激波冲击弦线中点（xi/a=0.5）的壁板颤振问题。并与Boyer 等[23］的数值结果进行对比。本研究中壁板在稳定颤振过程中壁板3/4 位置处振幅为0.459 6h，振动频率约为299.9 Hz。图2 为验证算例中壁面上的时均摩擦系数Cf分布情况，实线为本文数值计算结果，三角符号为Boyer 等[23］的计算结果。由于三维效应的存在，壁板3/4 处振动振幅和频率与文献[23］中三维计算结果（振幅0.472 1h，频率约294.84 Hz）存在一定的偏差，但壁面的时均摩擦系数仍能较好地吻合于文献结果，表明本文所采用的数值方法在激波冲击壁板计算中具有一定的准确性。

图2 壁板上的时均摩擦系数Fig.2 Time-averaged skin friction coefficient over flexible panel

2 计算物理模型

2.1 几何模型与边界条件

激波/边界层干扰作用下弹性壁板的流固耦合计算模型如图3 所示，壁板上表面受激波/边界层干扰作用，下表面为空腔，与文献[21，28］中的计算模型相同。壁板上表面入射斜激波前的马赫数、速度和密度分别为Ma1、U1和ρ1，弹性壁板密度、长度、厚度、振动位移、无量纲质量、弹性模量和泊松比分别为ρs、a、h、w、μs、E和v，其中0.3。无量纲动压定义为λ=ρ1U12a3D=875，抗弯刚度

弹性壁板上表面为激波/边界层干扰作用下流场计算获得的气动压力载荷，下表面空腔压力根据激波冲击位置计算得到，其计算式[21，28］为pc=(xi a)p1+[1-(xi a)]p3，其中xi为无黏条件下斜激波的冲击位置（壁板上任意一点离壁板前固定端的距离），p1和p2为入射斜激波前后气体压力，p3为反射斜激波后的气流压力，来流压力比p3/p1=1.4，对应入射激波角σ=35.58°。

来流和计算域上方部分边界设置为超声速入口边界条件，边界处的物理量给定来流参数，其中斜激波反射后的物理量采用激波关系式[36］计算。反射激波后的计算域上表面和出口均采用超声速出口边界，流动物理量直接由计算域内部外推得到。对于黏性流动，壁面采用无滑移绝热边界条件，基于壁板长度的来流雷诺数R ea=120 000，壁板前缘边界层厚度为δLE=0.015 6a。整个计算域长度为0 ＜x/a＜ 2.0，其中弹性壁板位于0.6 ＜x/a＜ 1.6。壁板两端设置为简支边界条件。此外，以激波/边界层干扰作用下刚性壁面的稳态流场结果作为流固耦合计算的初始值。

2.2 网格与时间步长无关性验证

计算区域采用结构网格划分，并在y方向对弹性壁板表面处的网格做加密处理。为保证计算结果的准确性，采用3 套网格（201×105、301×133、501×133）进行无关性验证，同时分析时间步长对壁板颤振特性的影响。网格和时间步长对弹性壁板3/4 位置处振动位移响应如图4 所示，其中时间t以a U1进行无量纲处理（下同）。由图4 可知，网格数量301×133、时间步长Δt=0.000 5a/U1工况下加密网格数量和减小时间步长对计算结果的影响较小，因此本文网格数量选择301×133，时间步长设置为Δt=0.000 5a/U1。

图4 网格和时间步对激波/边界层干扰下壁板3/4 位置处振动位移的影响（xi a=0.4）Fig.4 Effect of grid resolution and time step size on panel oscillation at3/4chord of flexible panel in shock/boundary layer interaction（xi a=0.4）

3 计算结果与分析

3.1 壁板颤振特性

与传统均匀超声速来流下壁板的气动弹性响应不同，激波/边界层干扰下壁板上表面承受较大的非均匀气动载荷，壁板上下表面压力差较大，弹性壁板将会偏离其初始平衡位置迅速开始振动。激波冲击位置xi a=0.4 工况下壁板的气动弹性响应如图5 所示。由图5（a）壁板3/4 处的位移时间可以看出，在初始压差的瞬态作用下，壁板逐渐偏离初始平衡位置，并开始变形振动。初始阶段振动位移较大，最大振动位移约为3.11h，随后振动位移逐渐衰减，经过若干瞬态振荡周期后，壁板振动位移持续增大并达到了稳定颤振状态。由图5（b）可见，壁板振动相图为封闭的环状曲线（其中ẇ为壁板振动速度），壁板的振动表现为极限环颤振状态。壁板不同位置处振动的频谱图如图6 所示，其中St=fa/U1为壁板颤振的无量纲频率，f为壁板颤振频率，壁板振动主频率St=0.299，同时选取壁板上不同位置的振动不影响壁板整体振动响应。由于壁板3/4 位置处的振动变形量较大，更易于监测，本文选取壁板3/4 位置处做为研究监测点。

图5 壁板3/4 位置处振动响应Fig.5 Vibration response at 3/4 chord of flexible panel

图6 壁板不同位置处振动频谱图Fig.6 Spectrum of vibration at different positions of panel

激波冲击位置xi a=0.4 工况下壁板颤振过程中的瞬态振形如图7（a）所示，其中Instant、Mean 分别表示壁板在一个周期内的瞬态值和时间平均值，Φ=0°，90°，180°分别表示x=1.35a处（壁板3/4 长度监测位置）振动位移的最大、中间和最小值，如图7（a）所示，壁板颤振呈现出近似2阶振动模态，壁板前后两端平均变形量相对于激波冲击位置并不对称，前端的振动幅值小于后端；此外，在一个颤振周期内，壁板总是围绕平均变形位置发生振动，壁板后端振动位移始终为负值，而壁板前端振动位移始终为正值。图7（b）为弹性壁板振动的x-t图，表明壁板振动呈周期性变化。

图7 壁板瞬态变形和平均变形分布与壁面振动x-t 图Fig.7 Mean and selected instantaneous panel shapes，and x-t diagram of panel deflections

壁板颤振Φ=180°时刻处的瞬时压力和密度云图如图8 所示，壁板前半部分主要发生“上凸”变形，因此其前端形成了前缘激波。而激波冲击位置下游壁板呈现“下凹”变形，壁板表面形成了明显的膨胀波和压缩波，边界层外围存在分离和再附着激波，同时边界层在入射激波冲击位置附近显著增厚。

图8 Φ=180°位置流场压力云图和密度云图Fig.8 Pressure contour and density contour for panel flutter at Φ=180°

3.2 颤振诱导的流动分离特性

流动分离是激波/边界层干扰的重要现象，通常发生在以激波主导的流动中。当入射激波足够强时，作用在流动上的逆压梯度会降低流动速度，边界层厚度增加，最终导致流动分离，并在激波入射点附近会生产一个分离区。xi a=0.4工况下壁板颤振对激波/边界层流动分离特性的影响如图9 所示。由图9（a）分离区长度随时间的变化过程可见，分离区长度在初始阶段变化较大，经过若干周期后，出现周期性振荡，其变化趋势和振荡频率与壁板振动过程保持一致，如图9（b）所示。当壁板达到稳定颤振状态时，分离区长度从Lsep/a= 0.288 振荡到Lsep/a= 0.376，即流动分离产生的分离泡随着壁板振动而发生周期性变化。以刚性壁板条件下，分离区长度对壁板颤振时的分离区长度进行无量纲处理（Lsep/LRigid），其中刚性壁板分离区长度约为0.324a，如图9（c）所示。当壁板处于Φ=180°位置，流场分离区长度缩减约11.1%；当壁板处于Φ=0°位置，流场分离区长度增加约16.05%。

图9 壁板颤振对流动分离的影响（xi a=0.4）Fig.9 Effect of panel flutter on flow separation（xi a=0.4）

刚性壁面工况和壁板颤振过程中的表面摩擦系数和压力分布如图10 所示。从图10（a）可以看出，壁板颤振过程中流动分离起始位置差异较小，基本位于x/a= 0.87 附近；而分离区的再附着点位置在一个周期变化较为明显，表明流动分离区域长度的变化主要是由再附着点位置变化引起。此外，与刚性壁面条件下的相比，流固耦合作用下壁板的颤振对激波冲击位置下游的压力分布和表面摩擦系数影响较为剧烈。

图10 壁板表面摩擦系数和瞬态压力分布Fig.10 Skin-friction along panel and instantaneous pressure distribution on panel

刚性和弹性壁板的分离区流场结构如图11所示。激波/边界层干扰下流动分离产生的回流区呈现狭长的扁平状，其流向尺度远大于法向尺度。壁板颤振对分离泡形状、再附着点位置、中心位置以及最大高度均有显著影响，在壁板颤振过程中，分离泡的形状被拉长，涡中心在激波入射点附近前后运动，在附着点向壁板后缘移动。通过对比图11 的图像，在Φ=180°时刻的分离区长度为0.376a，然而Φ=0°时刻的分离区长度仅有0.288a。文献[28］壁板强迫振动下激波/边界层干扰研究结果表明分离区变化主要与壁板表面曲率变化相关。针对本文xi a=0.4 工况，在激波冲击位置上游，壁板的表面曲率变化小，而激波冲击位置下游壁板表面曲率发生了显著变化，从而导致分离泡的形状和再附点位置发生了动态的变化。因此，可通过改变壁板的表面曲率，以便获得有利或不利的逆压力梯度，实现抑制或促进流动分离。

图11 工况xi a=0.4 刚性壁板和弹性壁板颤振的分离区流场结构（y 方向比例放大10 倍）Fig.11 Streamwise velocity contours for xi a=0.4 and over flexible panel and rigid-flat panel（Scale has been magnified by factor of 10 in ydirection）

不同工况下弹性壁板及刚性壁板不同轴向位置处的速度剖面如图12 所示，x/a=0.975，1.1位置处的轴向速度在0 ＜y a＜0.5×10-2范围内出现负值，这表明此处存在回流区域。由该组曲线可见，弹性壁板工况下，激波冲击位置上游区域（x/a＜1）的速度梯度高于刚性壁板工况，在附着点附近（x/a=1.225），弹性壁板工况下的速度梯度明显降低。

图12 沿计算域不同轴向位置的速度分布Fig.12 Velocity distribution along different axial positions in computational domain

3.3 激波冲击位置对壁板颤振特性的影响

不同激波冲击位置工况下弹性壁板3/4 处的位移时间响应历程如图13（a）所示，可以看出激波冲击位置对壁板颤振有显著影响；对于xi a=0.3 和xi a=0.7 工况，由于更靠近壁板固定端，壁板动力学响应经过若干时刻会收敛，呈现为静稳定状态，同时xi a=0.7 工况下壁板达到稳态所需的时间较短。xi a=0.3 工况相位图如图13（b）所示，可以清晰地看到壁板的动力学响应最终会收敛，这与xi a=0.4 时的极限颤振状态不同（图5（b））。

图13 激波冲击位置对壁板颤振的影响Fig.13 Effect of shock impingement location on panel flutter

图14（a）所示为激波冲击位置（xi a）对壁板3/4 位置处振动位移颤振特性的影响，xi a=0.5工况下振动位移幅值最大，当激波冲击位置逐渐偏离xi a=0.5 向壁板上游或下游移动时，壁板3/4 处的颤振振幅逐渐减小，而振动中心位置变化趋势并不单调，同时壁板颤振主频率逐渐减小。不同工况下壁板3/4 处振动位移的无量纲功率及功率谱密度如图14（b）和图14（c）所示，对于0.4 ≤xi a≤0.6 的工况，在St≈0.3 处均出现强烈的频率峰值，同时在St≈0.585 处出现较低的频率峰值。

图14 壁板振动功率谱分析Fig.14 Power spectrum analysis of panel flutter

初始时刻壁板上下表面压力差分布如图15所示，改变激波冲击位置，壁板所受到的综合压力载荷分布呈非线性变化，但压力差为0 的位置一直位于壁板的前端。同时，随着激波向下游移动，沿壁板前端的压力差分布逐渐平缓，而沿壁板后端的压力差分布逐渐剧烈。因此，壁板稳定颤振时前部分的变形量始终小于后部分。

图15 初始时刻壁板上下表面压力差分布Fig.15 Comparison of pressure difference between upper and lower surface of panel

激波冲击位置对振动平均变形形状的影响如图16 所示，不同冲击位置下壁板的平均变形总是呈现出2 阶模态，壁板后半部分的变形量始终大于前半部分，且变形量相对于激波冲击位置不对称。随着冲击位置偏离中心位置xi a=0.5，壁板前部分的振动平均变形量逐渐增加。但是，对于静稳定xi a=0.3 工况壁板后半部分的变形量最大，xi a=0.7 工况壁板前部分的变形量最小。

图16 壁板平均变形分布Fig.16 Mean deformation distribution of panel

为进一步研究壁板压力载荷与振动位移之间的关系，壁板3/4 位置处y方向的振动速度响应和壁面压力响应之间的关系如图17 所示。对于xi a=0.55，压力振荡与振动速度的相位基本一致，即壁板上表面受到脉动变化的压力载荷，壁板振动的加速度也呈现周期性变化，壁板振动速度和加速度同频发生改变，因此壁板处于稳定颤振状态。但是对于xi a=0.3，振动位移与压力响应之间存在相位差，这表明压力载荷对壁板振动具有反作用，导致系统振动能量衰减，最终达到静稳定状态。

图17 壁板3/4 位置处y 方向的振动速度和压力响应对比Fig.17 Comparison of vibration velocity in y-direction and pressure response at 3/4 chord of panel

3.4 壁板冲击位置对流动分离的影响

激波冲击位置对边界层流动分离区长度的影响如图18 所示。由图18（a）可见，随着激波冲击位置从xi a=0.35 开始向下游移动，分离区平均长度逐渐增加，分离区长度振幅变化较小。对于xi a=0.3 和xi a=0.7 工况，壁板振动最终呈现静稳定状态，由激波诱导的分离区长度同样处于稳定状态。图18（b）为激波冲击位置（xi a）分别位于0.35、0.4 和0.6 工况下分离区长度的时间历程，对于xi a=0.3 5 工况，与xi a=0.6 工况相比，分离区平均长度显著减小，振幅减小约26.65%，振荡达到稳定状态所需的时间也显著减小。同时，xi a=0.3 5 工况相较于xi a=0.6 工况，壁板3/4 处振幅减小约27.9%，导致激波/边界层干扰下分离区振荡得到抑制。

图18 激波冲击位置对流动分离的影响Fig.18 Effect of shock impingement location on flow separation

将不同激波冲击位置工况下弹性壁板的分离区长度由刚性壁板结果进行无量纲处理，结果如图19 所示，对于xi a＞0.4 的工况，流固耦合作用下激波/边界层干扰产生的分离区长度大于刚性壁板条件，即Lsep/LRigid＞1；对于xi/a=0.35 工况，壁板颤振过程中分离区长度始终小于刚性壁板条件的分离区长度，如图20 和图21所示，相对于刚性平板，壁板颤振过程中流动分离起始位置向下游移动，同时部分时刻再附着点位置向上游移动，二者共同作用下导致了分离区长度变短。因此，激波冲击位置位于xi a=0.35 时，壁板颤振能够有效抑制激波/边界层干扰流动分离。

图19 分离区长度相对变化值Fig.19 Separation length normalized by corresponding rigid-flat panel result

图21 xi a=0.35 工况下刚性壁板和弹性壁板颤振的分离区流场结构（y 方向比例放大10 倍）Fig.21 Streamwise velocity contours for xi a=0.35 and over flexible panel and rigid-flat panel（Scale has been magnified by factor of 10 in y-direction）

不同冲击位置下弹性壁板3/4 处Φ=180°时刻的压力云图如图22 所示，改变冲击位置对激波/边界层干扰的流动分离形态有着明显的影响。综上，随着激波冲击位置向下游移动，由激波诱导的分离泡的中心位置、再附着点位置也逐渐后移，这与图20 中的计算结果一致。在0.3 ≤xi a≤0.55 范围内，随着激波冲击位置后移，分离泡的形状被拉长、分离区最大高度显著增加。

图20 xi a=0.35 工况下壁面摩擦系数Fig.20 Skin-friction along panel with xi a=0.35

图22 不同工况下Φ=180°时刻的压力云图（y 方向比例放大10 倍）Fig.22 Pressure contours at Φ=180° under different working conditions（Scale has been magnified by factor of 10 in y-direction）

在Φ=180°时刻不同轴向位置速度沿y方向分布曲线如图23 所示。从x/a= 0.85 处的速度分布可以看出，随着激波冲击位置向上游移动，该位置处近壁面区域（0 ＜y a＜0.5×10-2）速度梯度减小，但对于xi a≥0.55 工况没有显著影响。通过观察壁板的速度分布，在0 ＜y a＜0.75×10-2范围内速度出现负值，这表明该位置出现回流区。对比壁板不同位置处的速度分布，不同激波冲击位置工况下处的回流区高度也受到影响。从x/a=1.35 处的速度分布可见，xi a=0.65 工况在该处的回流区高度最大，而xi a≤0.55 工况下的轴向速度（u u1）始终大于1，表明激波/边界干扰下流动分离在该位置处已结束。通过对比x/a=1.35 和x/a=1.7 处的速度分布发现，随着激波冲击位置向下游移动，壁板下游近壁面区域的速度梯度逐渐减小，但x/a=1.4 处的变化趋势较x/a=1.7 处更明显。