摘要
为探究民用飞机低速大迎角下平尾抖振流动特征和机理，基于应力混合涡流模拟方法，对典型民用飞机构型在流动分离状态下开展了高精度流动仿真研究。首先，以NACA 0015翼型大迎角分离流动仿真为例，验证了数值方法对分离流动的解析能力。其次，针对某型民机干净构型开展数值模拟，并与风洞试验数据进行对比。对比分析显示，数值计算得到的平尾平均压强系数和脉动压强系数在不同流动状态下均与风洞试验值吻合良好，且随迎角的变化趋势一致，验证了该方法对机翼分离流动干扰下的平尾流场仿真适用性。最后，基于仿真结果与试验数据的对比，深入分析了机翼分离流动对平尾气动干扰的机理及演化规律。分析结果表明：在抖振边界附近，机翼局部流动分离仅引起平尾较小幅度的抖振载荷；随着迎角进一步增大，强烈的机翼尾迹分离涡逐渐主导平尾流场，导致平尾抖振载荷显著增大。而当迎角继续增大，平尾逐渐脱离机翼尾迹影响区域，抖振载荷幅值增长放缓，此时平尾流场转由其自身的大迎角分离流动主导，并伴随流动主频向高频迁移。本研究揭示了机翼大迎角分离流动影响下平尾抖振载荷特征的演化规律，为民用飞机安全性能评估与设计提供指导。
关键词
Abstract
To investigate the horizontal tail buffet flow characteristics and mechanisms for civil aircraft at low speeds and high angles of attack (AOA), this study conducted high-precision flow simulations for a typical civil aircraft configuration under flow separation conditions based on the stress-blended eddy simulation (SBES) method. Firstly, the accuracy of the numerical simulation in resolving separated flow is validated using a NACA 0015 airfoil at high AOAs. Then, numerical simulations were performed for a clean configuration of a civil aircraft and compared with wind tunnel test data. The comparative analysis shows that the numerically predicted mean and fluctuating pressure coefficients on the horizontal tail agree well with the wind tunnel test data under various flow conditions and exhibit consistent variation trends with AOA, thus verifying the method's suitability for simulating the horizontal tail flow field under the interference of wing separation flow. Finally, based on the comparison between the simulation results and experimental data, the mechanism and evolution of aerodynamic interference on the horizontal tail caused by wing separation flow are analyzed in depth. Analysis results reveal that near the buffet boundary, local wing flow separation induces only minor buffet loads on the horizontal tail. As AOA increases, strong wing wake separation vortices gradually dominate the flow field of the horizontal tail, significantly amplifying the buffet loads of the horizontal tail. As the AOA increases further, the horizontal tail gradually moves out of the wing wake region, causing the growth of the buffet load amplitude to slow down. Consequently, the flow field becomes dominated by the horizontal tail's own high-AOA separation flow, accompanied by a shift of the dominant flow frequency toward higher frequencies. The study clarifies the evolution of buffet load characteristics of the horizontal tail under the influence of wing separation flow at high AOA, offering guidance for safety assessment and design of civil aircraft.
飞机在低速大迎角流动状态下，机翼尾迹区的分离流动会引起平尾抖振问题，严重影响飞机结构安全和操稳特性[1-3]。为研究飞机的抖振特性，准确地预测非定常气动载荷至关重要。目前，基于雷诺平均N-S方程（reynolds averaged Navier-Stokes，RANS）的数值仿真方法在预测附着流动的气动载荷方面表现良好，且对计算资源需求量较小；但对大范围分离流动问题，其预测精度会显著下降[4-5]。大涡模拟（large eddy simulation，LES）方法通过模拟小尺度涡，在大迎角分离流动问题解析精度较好，但LES方法对计算机资源需求巨大，这也限制了其在工程中的应用[6-8]。
目前，RANS-LES混合方法（如DES（detached eddy simulation）和SBES（stress-blended eddy simulation））结合了RANS高效率和LES高精度的特点，其在大规模工程湍流问题中得到广泛应用[4，9]。陈浩等[10]和刘周等[11]研究验证了DES类方法对NACA 0015翼型[12]大迎角失速分离模拟的精度和可靠性。肖志祥等[13]基于S-A湍流模型与LES方法构造了RANS-LES混合方法，证明了其对翼型段和椭球体大迎角绕流流动分离的模拟能力。宋科等[14]研究表明混合RANS-LES方法对多段翼型的最大升力系数与最大升力迎角具有较好的计算精度。白俊强等[15]通过混合方法对大迎角下的钝前缘大后掠三角翼模型进行了仿真研究，验证了DES和DDES（delayed detached eddy simulation）方法在解析复杂涡系演化方面的精度与能力。然而，针对飞机这类复杂外形的高雷诺数湍流问题，由于DES类方法需要更精细的网格才能解析LES区域，其计算成本仍比RANS方法模拟高出数个量级，这限制了其在多工况设计迭代和多学科耦合问题的应用。与DES类方法相比，近年来发展的应力混合涡流模拟（SBES）方法在RANS区域采用了更强的屏蔽函数[16]确保在近壁区域保持RANS模型的稳定性，并通过引入额外的源项，实现RANS和LES区域的快速转换，使其在相同网格分辨率下能够解析更丰富的湍流结构[17]，因而更适用于工程复杂外形的高精度湍流模拟。Ekman等[18]研究表明SBES方法对钝头体高雷诺数分离流动具有较好的计算精度，同时该方法对时间步长敏感度较低，能够在不损失精度情况下缩短仿真时间，提高计算效率。Ekman等[19]还评估了不同混合RANS-LES方法对钝头体不同工况下分离流的解析精度，验证了SBES方法在计算精度和捕捉流动物理细节方面具有显著优势。
综上所述，SBES方法对飞机复杂外形高雷诺数流动分离问题具有一定的适用性，但仿真精度还需通过风洞试验做进一步验证。本文应用SBES方法对飞机低速大迎角分离流动开展仿真研究，结合已有的风洞试验结果，对典型状态的飞机的脉动压强变化进行对比分析，研究SBES方法对抖振流动的仿真精度；并进一步研究不同状态抖振流动发展规律和流动机理。
1 数值方法与算例验证
本文中非定常仿真采用的应力混合涡流模拟方法[17]（SBES）属于混合RANS-LES方法的一种，其在近壁面采用RANS方程模型，而远壁面采用LES方程模型。该方法湍流应力张量的混合函数如下：
(1)
式中：、分别为雷诺应力张量项的RANS部分和LES部分，fSBES为屏蔽函数[16]。SBES方法还可以指定RANS模型和LES模型公式。如果RANS模型和LES模型基于涡黏假设，SBES方法的涡黏系数可由下式得到：
(2)
式中：、分别为湍流黏度的RANS部分和LES亚格子尺度模型部分。根据相关研究[13-14]，本文研究在RANS区域使用k-ω SST模型[20]，亚格子模型采用动态Smagorinsky模型[21]。数值仿真采用基于单元中心的有限体积求解器ANSYS Fluent 2022R1，其中动量对流项离散格式采用中心差分格式，其他对流项采用二阶迎风格式，梯度计算采用最小二乘算法。
针对NACA 0015翼型，弦长c=1.0 m，参考其相关风洞试验研究[12]，仿真条件为来流马赫数M∞=0.28，基于弦长的雷诺数Rec=1.95×106。计算网格采用非结构网格，翼型表面有407个单元点，并对尾迹区进行加密处理以保证对分离区的分辨率为1%c，壁面第1层网格高度为1×10-5c，保证壁面y+＜1。另外，三维网格为在二维网格的基础上沿展向拉伸0.5c得到[11]，以捕捉分离涡的三维效应，网格局部示意图见图1。翼型表面采用无滑移壁面条件，两端对称面采用周期边界条件，远场采用压力远场边界条件，流动采用全湍流计算。
为检验网格精度对仿真精度的影响，首先生成3套不同密度网格，分别在翼型表面和尾迹区进行加密，网格具体尺寸见表1。分别通过这3套网格对翼型分离流动进行仿真计算，时间步长固定为Δt=0.15 ms，对应的时间无量纲步长为δt=Δt·U/c=0.013（其中 U为来流速度），每个时间步子迭代步数需保证残差收敛到10-3数量级。仿真得到翼型α=16°状态下的平均升力系数 CL，mean，升力系数振荡幅度ΔCL以及升力系数振荡主导频率f的结果见表2。从表2可以看出，较稀疏网格A得到的平均升力系数和振荡幅值有较大偏差，较密的网格B和C平均升力系数和振荡幅度差别不大，不同网格振荡频率差异较小，因此认为网格B已满足仿真精度要求。
图1翼型计算网格
Fig.1Computational mesh of airfoil
表1不同网格尺寸参数
Tab.1Parameters of different mesh sizes
表2不同网格仿真结果
Tab.2Simulation results of different meshes
随后，通过网格B，采用不同时间步长仿真计算得到α=16°状态的非定常气动力结果见表3。从表3可以看出，时间步长δt=0.013的平均升力系数和升力系数振荡幅度仿真结果与时间步长δt=0.009的结果差别较小，因此可以认为时间步长δt=0.013已经达到仿真精度要求。
表3不同时间步长仿真结果
Tab.3Simulation results with different time steps
图2给出了RANS、SBES方法仿真计算得到的翼型段升力系数随迎角变化曲线，通过与文献[12]相同条件试验结果对比可以看出，RANS方法在小迎角计算与试验值吻合较好，在最大升力系数附近结果与试验值偏差较大。试验[12]测得翼型的最大升力系数出现在α=13°附近，而本文RANS方法得到的最大升力在α=15°附近，误差较大；应用SBES方法得到的最大升力也在α=13°左右，与试验值差距较小，且在失速后大迎角下的平均气动力变化趋势也更接近于风洞试验结果。
图2NACA 0015翼型升力系数随迎角变化（M∞=0.28）
Fig.2Variation of lift coefficient with AOA of NACA 0015 airfoil (M∞=0.28)
图3为α=16°状态下，不同数值方法计算得到的翼型表面平均压强系数分布与试验值[12]对比情况。从图3可以看出，两种方法计算的翼型上表面负压峰值均偏大，其中基于RANS方法计算结果偏大更为明显。在上表面负压区，SBES方法的计算结果比RANS方法更接近试验值。图4给出了两种数值方法得到的翼型α=16°状态下瞬时分离涡强度分布对比（这里使用Ω·c/U进行无量纲化，Ω为分离涡强度）。从图4可以得出，RANS方法未能捕捉到涡脱落现象，而SBES方法能够捕捉到涡脱落现象，这与刘周等[11]DDES方法仿真结果类似；另外，这也表明尾迹涡对下游流场的影响范围较广。综上所述，本文数值方法对翼型分离流的计算精度得到了验证，为后续飞机的抖振流动仿真奠定基础。
图3两种方法计算的翼型表面压强系数对比（α=16°）
Fig.3Comparison of pressure coefficients on airfoil surface from two methods (α=16°)
图4不同方法计算的大迎角下翼型流场涡量分布对比
Fig.4Comparison of vorticity distribution of airfoil flow field at high AOA computed by different methods
2 飞机大迎角分离流动仿真
2.1 仿真建模
本文针对某型民用客机风洞模型开展大迎角流动仿真分析，风洞模型为全机缩比模型（缩放比为5.5%），平均气动弦长c=0.44 m；模型在风洞中通过腹部单支撑机构[22]安装（图5（a）），上部分为模型和天平，下部分为迎角变化机构，可实现迎角变化范围为-8°~26°，精度为±3′。该风洞试验于中国航空工业空气动力研究院FL-9风洞中开展，在平尾上η=0.30~0.92区域不同站位布置了40个Kulite脉动压力传感器，其动态数据采样率为4 kHz，分布位置如图6所示，风洞试验具体细节详见文献[22]。该试验测得了不同飞机构型大迎角状态平尾的压强变化数据，用于平尾抖振载荷特性分析。
数值仿真采用风洞模型的干净构型（图5（b）），流动条件为M∞=0.3，Rec=3×106。计算网格采用非结构多面体网格，机翼和平尾表面网格分布如图7所示。根据翼型仿真建模研究，对机翼尾迹区域进行加密处理（图7（b）），保证机翼和平尾之间尾迹区域分离涡的捕捉精度满足Δx/c=1%，体网格单元数约2 000万。仿真计算在超算平台上开展，单节点CPU为128核心，内存为256 G。
基于SBES方法开展分离流动仿真计算，入口和出口分别为速度进口和压力出口边界条件，飞机表面为无滑移物面条件，对称面为对称边界条件。非定常时间步长为Δt=0.05 ms，对应的时间无量纲步长为δt=0.013，则气流流经一个平均气动弦长需要约77个时间步，每个时间步内迭代至少15次子迭代步数以保证残差至少下降到10-3数量级。计算结果达到稳定后，对飞机表面平均压强和脉动压强计算统计时间持续至少100个过流弦长时间（t·U/c），单个状态的仿真时间大约需要3 000个CPU核时。
图5风洞模型与仿真模型
Fig.5Wind tunnel model and simulation model
图6平尾上表面风洞试验压强测点位置分布
Fig.6Pressure gauge location distribution of wind tunnel test on the upper surface of horizontal tail
图7非定常仿真计算网格
Fig.7Computational mesh for unsteady simulation
2.2 平尾流场分布对比
图8给出了数值仿真得到α=15.0°~20.0°状态平尾上表面流场分布变化情况，包括平均压强系数Cp，mean和脉动压强系数Cp，rms分布。由图8（a）可以看出，当α=15.0°时，平尾表面负压区只分布在翼尖前缘较小范围；当α=17.5°时，平尾局部负压区向中部前缘附近移动，负压区的范围也有所扩大；当迎角增大到20.0°时，负压区向平尾根部和中间区域移动，且范围进一步扩大。由图8（b）可以看出，当α=15.0°时，脉动压强较大幅值（Cp，rms＞0.5）集中在平尾翼尖前缘区域。当α=17.5°时，脉动压强较大幅值向内侧移动，分布在平尾中部和外侧区域。当α=20.0°时，脉动压强较大幅值进一步向内侧中间区域移动，分布范围也进一步扩大。综上所述，平尾流场随迎角变化呈非线性规律，这与其所处的复杂流动环境密切相关。
随后，基于该飞机开展的低速大迎角风洞试验结果，对比分析本文仿真方法对平尾抖振流场的计算精度。图9为不同迎角下（15.0°，17.5°，20.0°），数值计算得到的平尾各个站位平均压强系数Cp，mean与风洞试验结果的对比，不同颜色点和线分别对应不同迎角下的试验值和计算值。由图9可以看出，当α=15.0°时，平尾η=0.30~0.60站位区域压强系数分布与试验差异很小，η=0.73~0.80站位前缘附近x/c=0.2较试验值偏小，中后缘与试验吻合良好。当α=17.5°时，不同站位计算值与试验值偏差较小。当α=20.0°时，平尾内侧η=0.50站位区域，前缘x/c=0.1处平均压强系数较试验值偏小；平尾外侧η＞0.50站位区域，前缘压强系数较试验值略微偏大，弦向变化趋势与试验一致。总体上看，随着迎角的增大，平尾内侧前缘负压峰值逐渐增大，尤其是平尾中间区域变化较为剧烈；计算值与试验值都反映了平尾外侧前缘负压峰值逐渐减小的变化规律。
图10为不同迎角下，数值仿真得到的平尾各站位的脉动压强系数Cp，rms与风洞试验结果的对比情况。
图8平尾表面流场分布随迎角变化
Fig.8Variation of flow field distribution of horizontal tail surface with AOA
图9不同迎角下，平尾各个站位的平均压强系数Cp，mean与试验值对比
Fig.9Comparison of mean pressure coefficient Cp, mean at various positions of horizontal tail with different AOAs with test data
图10不同迎角下，平尾各个站位的脉动压强系数Cp，rms对比
Fig.10Comparison of fluctuating pressure coefficient Cp, rms at various positions of horizontal tail with different AOAs
由图10可以看出，当α=15.0°时，平尾中内侧站位（η=0.30~0.60）只在前缘区域（x/c＜0.1）有较大的压强脉动，弦向中部和后缘区域脉动压强较小（Cp，rms＜0.1），这与试验值表现一致；平尾外侧站位（η=0.70~0.80）有较大范围的压强脉动，这区域计算值略小于试验值，在弦向上变化趋势与试验一致。当α=17.5°时，计算值与试验值较为接近，平尾中间站位η=0.40~0.60区域压强脉动范围明显增大，外侧站位η=0.73~0.80区域峰值增大同时，范围扩展到靠后缘附近。当α=20.0°时，内侧区域脉动压强峰值进一步增大；外侧区域脉动压强峰值减小，但后缘附近脉动压强增大，呈现出随弦向位置变化较小的趋势。当α=20.0°时，翼尖站位脉动压强计算与试验吻合良好，在平尾内侧η=0.30~0.40脉动压强峰值较试验值偏小，中间区域η=0.51~0.60，脉动压强仿真计算值较试验值有一定偏差，但总体上脉动压强系数随弦向位置变化规律与试验值保持一致。
综上所述，本文仿真方法精确地解析了机翼尾迹影响下平尾平均压强和脉动压强的分布变化规律，因而能够满足平尾抖振载荷精细化分析需求。
2.3 平尾抖振载荷特性分析
图11为仿真得到飞机平尾不同状态下升力（载荷）系数振荡幅度变化情况。
图11平尾升力系数振荡幅度变化（α=10.0°~20.0°）
Fig.11Variation of oscillation amplitude in lift coefficient of horizontal tail (α=10.0°–20.0°)
由图11可以看出，α≤12.5°状态下，平尾抖振载荷幅度ΔCL处于较低水平，说明此状态下机翼尾迹对下游的平尾影响较小，风洞试验结果也表明此状态平尾测点无明显压强脉动现象。当迎角增大到15.0°时，抖振载荷幅度ΔCL逐渐增大到12.5°状态抖振载荷的4倍，主要是由于此状态机翼尾迹区对平尾影响较大，前文平尾流场分布对比中的脉动压强分布也印证了这一点。当迎角增大到17.5°时，平尾抖振载荷幅度进一步增大，说明尾迹区对平尾影响继续增强。而当迎角达到20.0°时，抖振载荷有减小的趋势；这也与平尾流场分布对比中的脉动压强分布变化一致，平尾外侧区域（η＞0.60）在20.0°迎角时，脉动压强明显减小，说明此状态尾迹区对平尾外侧影响有所减小，因而抖振载荷幅度有所降低。
进一步分析平尾抖振载荷频率特征，图12给出了两种迎角（α=15.0°，α=20.0°）状态下，数值计算得到的飞机机翼和平尾的气动力系数随时间变化情况。由图12可以看出，随着迎角的增大，机翼平均升力系数减小，主要原因是迎角增大使得飞机分离区扩大，造成机翼失速效应增大。平尾平均升力系数随迎角增大而增大，这是由于机翼下洗流动影响，平尾较机翼更晚进入失速状态；另外，平尾升力振荡幅度也随迎角增大，这与机翼尾迹强度和平尾自身流动分离区域有关。
图12机翼和平尾升力系数时间历程
Fig.12Lift coefficient time history of wing and horizontal tail
随后，通过快速傅里叶变换方法得到两种迎角下机翼和平尾的升力系数频谱（图13）。由图13可以看出，当α=15.0°时，机翼主频峰值在60 Hz附近，在110 Hz附近也有幅值较小的峰值；平尾在不同频率有峰值，其中60、80、110 Hz附近频率的幅值接近，说明机翼在110 Hz的非定常流动对平尾的影响较大。当α=20.0°时，机翼上在40、80、110 Hz附近幅值接近，而平尾流动在110 Hz频率附近占据主导。
图13机翼和平尾升力系数频谱
Fig.13Lift coefficient frequency spectra of wing and horizontal tail
通过风洞试验数据做进一步分析，图14为不同状态下机翼后缘和平尾前缘测压点的压强变化试验值的频谱图。小迎角下（α=15.0°），机翼测点压强变化主频集中在80 Hz，而平尾测点压强变化在80、110 Hz附近幅值接近。大迎角下（α=20.0°），机翼测点压强变化在40、60 Hz幅值较大，在80 Hz幅值较小；平尾测点压强变化在110 Hz附近有较大幅值，在40~60 Hz范围幅值较小。综上所述，随着迎角的增大，机翼后缘压强的峰值频率向较低频区域移动，而平尾前缘压强的峰值频率向较高频区域移动。这种变化规律与仿真计算得到平尾升力变化规律一致，也进一步说明仿真方法准确捕捉到不同流动状态下平尾抖振频率特征发展规律。
为探究平尾抖振载荷特征随流动状态变化的流动机理，图15给出了不同迎角状态下仿真得到平尾内、外侧两截面的瞬时分离涡强度分布图，其中两个截面分别为平尾内侧截面（η=0.51）和外侧截面（η=0.73）。
由图15可以看出，α=12.5°时，机翼内侧无明显分离流动，尾迹区对平尾的影响较小；机翼外侧开始出现前缘分离，尾迹区对平尾有一定的影响，但此时机翼后缘涡强度较小且在机翼下洗流动影响下平尾当地迎角较小，因而抖振载荷处于较低水平。当时，机翼内侧也开始出现较大流动分离，使得平尾内侧脉动压强有明显增大；机翼外侧后缘分离涡强度进一步增大，对平尾的影响增强，这使得脉动压强幅度增大，引起抖振载荷幅度明显增大。随着迎角增大到α=15.0°时，虽然机翼后缘分离涡强度增大，但此时平尾外侧已经退出机翼尾迹和机翼下洗的影响区域；此时平尾当地迎角约为来流迎角，因此平尾自身大迎角状态导致了其分离流动。平尾自身分离流动主导的抖振载荷比尾迹区影响主导的抖振载荷小，同时平尾自身分离流动频率也由较高频率主导；这也解释了前文中平尾抖振载荷幅度和频率特征随迎角状态的变化规律。
图14飞机不同测点压强变化风洞试验值频谱
Fig.14Frequency spectra of pressure variation at different measurement points on the aircraft from wind tunnel test
图15飞机平尾内、外侧截面不同迎角下瞬时分离涡强度分布
Fig.15Distribution of instantaneous separation vortex intensity at inboard and outboard cross-sections of aircraft's horizontal tail at different AOAs
3 结论
本文通过基于应力混合涡流模拟（SBES）方法的高精度湍流仿真结合风洞试验结果对比验证，揭示了民机机翼分离流动影响下平尾抖振载荷特征发展规律和流动机理，得出结论如下：
1）平尾平均/脉动压强系数的仿真结果与风洞试验结果吻合较好，且迎角变化趋势一致性显著，验证了SBES方法能够较为准确地捕捉机翼复杂流动影响下的平尾抖振载荷特征。
2）不同迎角状态，机翼尾迹对平尾的影响强度不同：在抖振边界附近，机翼局部流动分离引起的抖振载荷幅度较小；随着迎角的增大，较强的机翼尾迹涡主导平尾流场，导致抖振载荷呈显著增长趋势。
3）飞机在较小迎角时，平尾流场由机翼尾迹涡影响主导，流动主频集中于较低频率；而在较大迎角下，平尾局部区域逐渐退出机翼尾迹涡影响范围，因此减缓了抖振载荷的增长，此时平尾流场转由其自身的大迎角分离流动主导，并伴随流动主频向高频迁移。
