张馨元
摘 要：采用Roe通量微分分裂格式离散Navier-Stokes（N-S）方程组的对流项，用统一时间步长的伪时间子迭代隐式LU-SGS方法进行时间推进，利用变形网格方法求解任意拉格朗日-欧拉（ALE）坐标下的N-S方程组，对NACA-0012翼型小幅振荡和动态失速以及三段翼型振荡襟翼等非定常流动进行了数值分析；通过与可获得的实验结果对比，吻合良好。与此同时，还对比分析了代数B-L、两方程SST及k-g等三种湍流模式的表现，总体上SST和k-g模式结果优于B-L模式；在动态失速流动中，SST模式较k-g模式更靠近测量值；采用SST模式对三段翼型振荡襟翼流场进行详细研究，分析流动特点及其与升力特性之间的关系。
关键词：动网格 深失速 振荡翼型/襟翼
中图分类号：V21 文献标识码：A 文章编号：1674-098X（2017）10（b）-0014-06
Abstract：The Roes Flux Difference Splitting is employed as the spatial discretization scheme for the convective terms of the Navier-Stokes equations. The fully implicit lower-upper symmetric- Gauss-Seidel with pseudo time sub-iteration is applied as the temporal marching methods for both N-S equations and turbulence model equations. Flows around pitching airfoils or flap are investigated by solving N-S equations based on a moving grid system which described by Arbitrary Largangian-Eulerian method. Three turbulence models， such as algebraic Baldwin- Lomax， two-equation shear stress transport （SST） and k-g models， are employed to study the flows around pitching and dynamic stall NACA-0012， and the pitching flap of the three-element airfoil. The SST and k-g models perform better than B-L model. In the dynamic stall case， the SST model is better than the k-g model. Flowfields and flow characteristics of the pitching flap are investigated in detail by SST with respect to lift coefficients.
Key Words：Moving grid； Dynamic stall； Pitching airfoil/flap
随着计算机技术和计算方法的飞速发展，计算流体力学（CFD）在飞机和导弹设计中得到了越来越广泛的应用。具有运动边界的非定常复杂流场，如机翼颤振、操纵面嗡鸣、绕流片摆动、襟翼偏转、外挂物投放等，一直是计算流体力学研究的热点和难点问题。若上述流动发生在跨音速阶段，则还需考虑激波/边界层的相互干扰；若飞机经历动态失速，需要准确模拟大范围的流动分离；当飞机起飞和着陆时，襟翼偏转或振动对飞机的升力、阻力以及俯仰力矩特性产生重要影响，若严重时可影响飞行安全和误导驾驶员操纵，有必要深入研究俯仰翼型/襟翼流动。
当数值求解Navier-Stokes方程组时，需要引入湍流模式使方程组进行封闭。非定常动边界问题时，尤其是激波/边界层干扰及大范围流动分离问题，湍流模式的作用更加明显，其表现不容忽视。在飞机气动性能分析中，代数B-L、一方程S-A模式及两方程SST模式的应用极为广泛；然而，在分析非定常动边界绕流时，很多研究仍只考慮B-L模式。
本文发展了自行研究开发的多块结构网格程序NSAWET[1]（N-S Analysis based on Window- Embedment Technology），求解随网格点一起移动的运动坐标系即Arbitrary Lagrangian-Eulerian （ALE）的N-S方程组作为控制方程，利用变形网格方法对振荡翼型/襟翼进行数值研究；采用3种湍流模式对NACA0012小幅振荡和动态失速，以及三段翼型振荡襟翼等非定常流动进行数值研究。
1 数值方法
基于动网格技术研究俯仰翼型/襟翼的非定常流动特点，因此变形或运动边界对控制方程组、湍流模式、时间推进方法以及边界条件的影响被作为本文工作的重点，现分述如下。
1.1 控制方程组
本文所采用的控制方程组为随网格点一起移动的运动坐标系即ALE的N-S方程组[2，3]，与Euler描述下的N-S方程组相比，ALE的N-S方程组对流项中增加了由网格运动而引起的通量。
对流项采用通量微分分裂类Roe格式离散，粘性项采用二阶中心差分格式离散。
1.2 湍流模式
主要采用B-L、两方程k-ωSST以及k-g模式使N-S方程组封闭。代数B-L模式形式简单，勿需求解额外方程，SST及k-g模式需要求解关于k、ω或g的湍流模式方程组。endprint
2 计算结果及分析
2.1 NACA0012翼型小幅振荡
振荡NACA0012翼型具有较为详细的实验数据，被许多研究者所引用，也被作为本文的验证算例。翼型绕1/4弦点的振荡方式为：α（t）=αm+α0sin（2kt），其中αm为平均攻角，α0最大振幅攻角，k=为无量纲的振荡频率，为振荡频率，c为弦长，为来流速度。采用C型网格，翼型弦长为1，周向和法向网格数目分别为501和101，法向第一层网格距离物面为1.0×10-5。
2.1.1 亚音速情形
根据文献[6]，振荡翼型的流动参数为：Ma=0.60，Re=4.8×106，αm=4.86°，α0=2.44°，k=0.0810。无量纲时间步长0.05，一个振荡周期约800步，内迭代20次，湍流模式为SST。
图1是升力系数和俯仰力矩系数随攻角变化的迟滞曲线及其与实验数据以及文献[6]的结果对比，图2是几个典型时刻对应的压力系数Cp曲线。从图1与图2可以看出，本文数值计算结果与实验吻合非常好，且升力系数比文献[6]中利用k-ε模式的结果更接近实验值。通过Cp分布可以看出，振荡过程中，上扬和下俯时激波的表现差异非常大，即便下俯时的瞬间攻角较小（α=3.49°），激波仍很强，可是翼型上扬时瞬间攻角较大（α=4.28°）激波反而较弱，从客观上反映了流动的迟滞效应。当流动中出现激波并逐渐增强后，激波与边界层的干扰非常严重，导致了流动的强烈非线性，增加了数值分析的难度。
2.1.2 跨音速算例
该算例[3，5，6]中流动参数为Ma=0.755，Re=5.5×106，αm=0.016°，α0=2.51°，k=0.0814。
图3给出了B-L、SST及k-g模式计算的升力和力矩系数随攻角变化的迟滞曲线，及其与实验及文献[6]参考值的对比。可以看出，无论是文献还是本文方法，所得到的升力系数均比实验偏低。从图3还可以看出，SST和k-g模式对升力和力矩系数的模拟均较好，而且相互之间的差别非常小。B-L模式所得的升力迟滞环稍小，而力矩迟滞环稍大，计算效果较差。
2.2 NACA0012翼型动态失速
动态失速是指运动的压力面在超过其临界迎角时流场发生非定常分离和失速的现象。根据粘性作用的大小，动态失速可分为轻失速和深失速。轻失速除表现出一般静态失速的特征外，还有非定常分离的强粘性/无粘相互作用性质，边界层厚度为翼型厚度量级。深失速则表现为全粘性现象，存在高度非线性的压力脉动和翼型表面上有大尺度旋涡运动，边界层厚度可达到弦长量级[7]。
采用两方程SST和k-g模式分别计算深失速问题。文献[7]给出了NACA0012翼型正弦振荡深失速的实验结果，振荡规律与2.1节所述相同，只是参数有所改变，Ma=0.20，Re=1.0×106，αm=15°，α0=10°，k=0.15。为准确分析流动中的非定常效应，更细致的描述流场，无量纲时间步长为0.005，因此每一时间步攻角变化相对较小。
升力和力矩系数随攻角变化的迟滞曲线见图4，在翼型上扬过程中，计算结果与实验符合较好；当翼型下俯运动时，流动中由于受到强逆压梯度的影响，翼型上表面分离剧烈，两种模式均不能很好的模拟流动的分离特点，对其非线性的力和力矩变化值模拟能力较差，然而可喜的是，他们均可不同程度地模拟下俯过程中力和力矩的变化趋势。从图中可以看出，SST比k-g模式效果更好，这得益于SST模式的弱非线性涡粘系数的定义和对逆压梯度流动中涡粘系数进行限制，使其具有剪切应力输运的特点。B-L模式不适于计算此类流动，因而没有采用。
2.3 三段翼型振荡襟翼
选择文献[8]中的三段翼型，其前缘缝翼和后缘襟翼偏角均为30°。采用法向外推的单块O型网格（如图5所示），网格点数为582×131，法向第一层网格到物面的距离为1.0×10-5。
首先，对马赫数0.2，攻角8.109°，Re=9.0×106的定常绕流采用SST模式进行数值分析，三段翼型表面的压强系数Cp对比如图5右所示，计算结果与实验数据[8]吻合很好。
其次，在上述定常流动的基础上，数值分析振荡襟翼的流动特点。给定襟翼绕其10%弦点振荡的运动方式为α（t）=αm+α0sin（2kt），其中k=0.10，αm为襟翼基本偏角（30°），α0为10°。计算得到的三段翼型的总升力系數与各部件升力系数随襟翼偏角变化迟滞曲线如图6所示；与之对应，图7给出一个周期内典型时刻的流线。对照图6与图7可以看出，一个周期内流动状态可以根据襟翼偏角分为以下几个阶段：
（1）20.0°到36.3°（图中8→1→2），总升力系数Cl与襟翼偏角基本保持线性变化，襟翼的升力系数呈弱非线性增加，斜率逐步减小；偏角较小时流动保持附着，当偏角较大时，在襟翼后缘逐渐出现小的分离泡；襟翼提供的升力在36.3°到达最大。
（2）36.3°到38.2°（图中2→3），总升力继续增大，直至到达最大升力点，但上升斜率下降很快；主翼和前缘缝翼升力继续增加，由于襟翼后缘出现了较明显的旋涡，襟翼提供的升力已有所下降。
（3）38.2°到40.0°（图中3→4），总升力随偏角的增加不升反降；襟翼升力急剧下降，后缘分离区持续增大；主翼和前缘缝翼的升力均经历了从升高到低的过程。
（4）40.0°到33.2°（图中4→5），由于襟翼偏角减小，三段翼型的各段升力均降低，襟翼后缘分离区达到最大，襟翼提供的升力达到最小；但是由于襟翼的俯仰运动所造成的扰动传播尚未返回，因而主翼和前缘缝翼的升力系数还处于下探过程中。
（5）33.2°到20.0°（图中5→8），襟翼升力随偏角减小反而升高，然而受到襟翼俯仰运动的影响，总升力减小；襟翼后缘分离区逐渐减小，直至重新达到附着状态。endprint
由上述分析可知，各部件升力系数变化均具有较强的非线性特性，其最大值不发生在最大偏角，襟翼升力最小值也不在最小偏角。虽然主翼对在升力的贡献占了主要部分，但非线性特性主要由襟翼运动引起分离区变化所致。前缘缝翼升力则一直保持较弱的非线性特征，由于和襟翼相隔较远，对襟翼的俯仰变化敏感程度低。
3 结语
基于变形结构网格，采用三种不同表现的湍流模式对振荡翼型/襟翼流动进行了数值分析。通过与标准算例的验证，表明本文方法的正确性和计算的可靠性；通过计算结果发现，两方程SST和k-g模式，尤其SST模式，能捕捉流动中的非定常效应和模拟较强逆压梯度，在分析运动边界流动时获得了与实验吻合一致的效果。通过对单段翼型小幅度、深失速等俯仰振荡的准确模拟，并将其扩展到三段翼型绕流中的襟翼俯仰运动，可获得增升装置中单独翼对整个机翼气动特性的认识和了解，为设计增升装置打下坚实的基础。
参考文献
[1]Chen Haixin， Zhang Yufei， Xiao Zhixiang， Fu Song， NSAWET Code and Its Application [A].5th Asia Workshop on Computational Fluid Dynamics[C].2006.
[2]I.Demirdzic， M. Peric， Space conservation law in finite volume calculations fluid flow [J].International journal for numerical methods in fluids，1988（8）：1037-1050.
[3]C.Michler，H.De Sterck，H.Deconinck，An arbitrary Lagrangian Eulerian formulation for residual distribution schemes on moving grids [J].Computers & Fluids，2003（32）：59-71.
[4]肖志祥.復杂流动Navier-Stokes方程数值模拟及湍流模型应用研究[D].西安：西北工业大学，2003.
[5]A.Jahangirian， M. Hadidoolabi， Unstructured moving grids for implicit calculation of unsteady compressible viscous flows [J].INTERNATIonAL JOURNAL FOR NUMERICAL METHODS IN FLUIDS，2005（47）： 1107-1113.
[6]Zhang L. P.，Wang Z. J.，A Block LU-SGS implicit dual time-stepping algorithm for hybrid dynamic meshes[J].Computers & Fluids， 2004（33）：891-916.
[7]钱炜祺，符松，蔡金狮.翼型动态失速的数值研究[J].空气动力学学报，2001，19（4）：427-433.
[8]Daryl L. Bonhaus， W. Kyle Anderson， Dimtri J. Mavriplis， Numerical study to assess sulfur Hexafluoride as a medium for testing multielement airfoils NASA technical[Z]. 1995.endprint
