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

    基于流固耦合的激波/边界层干扰作用下壁板颤振特性    刘为佳，李映坤，陈雄，李春雷南京理工大学 机械工程学院，南京 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 为控制体表面法向量，以上各式的具体形式Reference[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 panel2 计算物理模型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）可见，壁板振动相图为封闭的环状曲线（其中ẇ为壁板振动速度），壁板的振动表现为极限环颤振状态壁板不同位置处振动的频谱图如图6 所示，其中St=fa/U1为壁板颤振的无量纲频率，f为壁板颤振频率，壁板振动主频率St=0.299，同时选取壁板上不同位置的振动不影响壁板整体振动响应由于壁板3/4 位置处的振动变形量较大，更易于监测，。