在跨声速机翼表面存在着明显的附面层分离和激波与附面层干扰等复杂的流动现象。跨声速状态下激波附面层干扰引起的跨声速抖振问题具有复杂的非定常、非线性特征，极大地制约了跨声速飞行器的飞行包线[1]。抖振机理的研究有利于揭示激波自激振荡的机制、激波附面层干扰的机理以及压力脉动的产生等，从而实现对激波抖振的有效控制。1976年，McDevitt等[2]通过实验系统研究了激波抖振现象，发现对于厚度18%的双圆弧翼型，在一定的马赫数下，激波会沿着翼型上下表面周期性往复运动。Tijdeman和Seebass[3]总结了以往的跨声速激波抖振的研究，将其分为A、B和C三类，其中绕厚度18%对称双圆弧翼型的流动即为C型，激波在翼型上下面传播，向上游运动的过程中逐渐减弱，并最终离开翼型前缘。2001年，Lee[4]提出了一种自激反馈模型，对激波抖振的机理进行解释，认为激波的周期运动诱导的压力波在附面层的分离区中向下游传播，到达尾缘时产生的声波又向上游传播，为激波的振荡提供能量，维持系统的稳定。Lee的模型得到了许多的研究者[5-10]的验证，对以后的研究起着重要的指导作用，但是其对后缘处声波的产生未给出确切解释，仍需进一步研究。
模态分析方法能够对复杂高维流场进行低维分解，捕捉对于非定常流动发展起重要作用的模态，并在时间和空间尺度上研究流动的发展。近年来，基于流场特征提取的模态分析方法，本征正交分解(POD)[11]和动态模态分解(DMD)[12]，常被应用于复杂非定常流场的分析中。潘翀等[13]运用DMD方法从时间和空间维度上分析了带襟翼翼型尾流的流动特征，捕获了尾流场基波和高阶谐波的频率、波长和传播速度等丰富的信息。刘宏康等[14]详细研究了开口空腔的非定常流动，结合POD和DMD方法分析了流场的压力脉动特征，揭示了空腔后缘声波辐射的2种不同路径。寇家庆等[15]基于POD和DMD方法，研究了A型跨声速激波抖振的压力场，对比了2种模态分解方法的特点。陈立为等[10]则研究了双圆弧翼型激波的周期性运动，提出了针对C型激波抖振的自激反馈模型，并用POD分析了主要模态与激波运动以及尾缘压力脉动的关系。但基于系统能量排序的POD方法并不能准确捕捉与抖振频率相关的主要模态，而基于频率排序的DMD方法在保留POD方法优势的同时，更包含了系统的动力学特性，能提取流场的各阶特征频率。本文针对厚度18%的对称双圆弧翼型，在时间维度上对于跨声速瞬时流场运用DMD方法，提取包含流场各阶主频的模态，分析了流场的动力学特性。
1 数值方法本文计算使用基于SA(Spalart-Allmaras)一方程湍流模型[16]的延迟分离涡模拟(Delayed Detached-Eddy Simulation, DDES)[17]求解绕厚度18%对称双圆弧翼型的非定常流动，重构格式为三阶MUSCL(Monotonic Up-stream-centered Scheme for Conservation Laws)[18]，并使用隐式双时间步推进[19-20]。
2 DMD方法DMD基于流场的动力学特性，提取的各阶模态在时间和空间演化特性上相互正交。本文截取流场进入稳定发展的非定常状态后的500个样本进行模态分解，采样频率为5 555.6 Hz，后续讨论样本数量的关系，也基于此样本。
取瞬时流场矩阵U1N-1和U2N为
|
(1) |
式中：ui为第i个瞬时流场数据列向量；N为所取的瞬时流场样本总数目。相邻瞬时流场之间的时间间隔为Δt，假定Δt很小，则两相邻瞬时流场满足线性映射关系：
|
(2) |
通过求解矩阵A的特征值和特征向量，以求得相应流场的特征信息。实际计算中，一般求解矩阵A的低维近似矩阵F代替求解高维矩阵。求解方法为对U1N-1进行奇异值分解，即
|
(3) |
式中：W和V均为酉矩阵；Σ为对角矩阵。再求近似矩阵F，使得A=WFWH，近似矩阵F可通过求解式(4)得到：
|
(4) |
即使得矩阵U2N-WFΣVH的Frobenius范数最小。从而得到：
|
(5) |
对矩阵F进行奇异值分解，得到其特征值μj和特征向量Λj，μj即为Ritz特征值[21]，其对数形式的实部和虚部分别表示增长率gj=Re(lg μj)/Δt和频率wj=Im(lg μj)/Δt。对应的0动态模态Φj和模态能量||Φj||[12]为
|
(6) |
图 1给出了厚度18%对称双圆弧翼型物理模型及网格，翼型弦长c=203 mm，展向长度为0.2c，远场为40c。本文使用多块对接结构网格，壁面设置无滑移绝热壁，展向为周期边界，远场为压力远场边界。使用两套网格计算，grid1、grid2网格分辨率周向、法向和展向分别为465×171×41、625×201×41，第1层网格高度使得y+ < 1。
|
| 图 1 对称双圆弧翼型壁面及对称面网格 Fig. 1 Symmetric circular-arc airfoil wall and symmetry plane mesh |
本文计算来流马赫数Ma=0.76，基于弦长的雷诺数Rec=1.1×107，迎角为0°，无侧滑，物理时间步长为0.002c/a∞，a∞为来流声速。图 2给出了翼型上表面时均压力系数Cp沿流向分布，grid2的计算结果与文献[2]和文献[10]结果基本吻合。图 3给出了升力系数的功率谱密度(PSD)曲线，使用最大熵功率谱估计方法，在幅值上本文计算结果稍大，但变化趋势与文献一致，主频相近，本文计算得到的斯特劳哈尔数St为0.150(f=StU∞/c, f为脉动频率，U∞为来流速度)，文献[10]结果为St=0.148，该主频对应激波的自激振荡频率。
|
| 图 2 壁面压力系数分布 Fig. 2 Wall pressure coefficient distribution |
|
| 图 3 升力系数功率谱密度 Fig. 3 Lift coefficient power spectral density |
对所取样本进行DMD分析，得到流场的各阶主模态。图 4给出了样本数N=100、200、500下的DMD谱，对应的采样频率分别为1 388.9、2 777.8、5 555.6 Hz。随着采样频率升高，捕捉流场特征模态的频率也越高，且样本数200和500的结果趋近，本文取样本数500的结果进行分析。图 4中空心方框为本文所取的前5阶模态，其中频率为0的模态为静模态，表征流场的平均特性，其余模态均成对出现，其特征值为共轭复数。所取模态增长率/衰减率均在0附近，说明对流场的非定常发展过程具有持续作用。
|
| 图 4 频率与增长率/衰减率的关系 Fig. 4 Relationship of frequency with growth rate/decay rate |
图 5给出了DMD的Ritz值以及各阶模态能量与频率关系。各阶模态的特征值均处在单位圆附近，个别模态位于单位圆内，说明所取样本处于准中性稳定状态。以往研究者通常根据模态幅值或者模态能量进行排列和选取特定的模态。本文根据后者选取图 5中前5阶模态进行分析，各阶模态主频分别为激波自激振荡频率的倍数，对激波的自激振荡起主要作用。
|
| 图 5 DMD的Ritz值和模态能量与频率关系 Fig. 5 Relationship of Ritz value and mode energy with frequency of DMD |
图 6给出了DMD前4阶模态(Mode 1~Mode 4)的模态系数a(t)随时间t的变化及其功率谱密度曲线。结果表明，前4阶模态系数具有简谐振动特征，Mode 1幅值随时间衰减，Mode 2~Mode 4幅值则随时间增长，与图 4中的模态的增长率/衰减率一致，同时模态系数的幅值随阶数的升高而减小。对各阶模态系数进行傅里叶分析，各阶模态包含单一主频，其对应的St分别为0.141、0.295、0.436、0.577，对应激波自激振荡St的约1~4倍，且幅值逐渐减小，其中Mode 1在激波抖振过程中占主导作用。
|
| 图 6 DMD前4阶模态系数随时间的变化及其功率谱密度曲线 Fig. 6 Variation of coefficient of the first four modes of DMD with time and its power spectral density curves |
图 7(a)~(d)分别给出了Mode0~Mode3的压力脉动(实部)的空间分布，p∞为来流压力，文献[22]研究表明实部和虚部在流动特征上区别不大, 图中无量纲压力值的正负代表着正负压力脉动。DMD各阶主模态阶数越高，空间维数也越高，在翼型上下面，正负压力脉动相间分布，奇数阶主模态压力脉动反对称分布(上下面压力脉动符号相反幅值相等)，偶数阶主模态压力脉动对称分布(上下面压力脉动符号相同幅值相等)。压力脉动的区间主要集中于翼型的0.5c~0.85c处，即激波抖振的区间，与文献[10]结果相吻合。
|
| 图 7 DMD模态的空间分布 Fig. 7 Spatial distribution of DMD modes |
圆弧翼型的激波自激振荡属于C型激波抖振，激波在翼型上下面周期性往复运动，激发压力波的传播。取t0、t0+1/6T、t0+2/6T、t0+3/6T 4个时刻观察激波抖振的过程，其中t0为激波自激振荡某个周期的初始时刻，T为振荡周期。图 8给出了激波运动半周期内对称面数值纹影图和DMD第1阶模态随时间的演化过程，Δρ为密度梯度。观察图 8纹影图激波的运动，对称翼型，激波在翼型上下面运动的过程是一致的，单独分析其在翼型上面的运动, 半个周期内，激波从靠近翼型尾缘开始向上游运动至翼型中部，运动区间大致为0.5c~0.85c，如图 8纹影图所示。初始时刻激波分叉较多，激波较强，如图 8(a)纹影图所示；向上游运动过程中，激波逐渐变弱，激波后分离区逐渐增大，至翼型中部时，激波压缩成一道，分离区扩大至最大，如图 8(d)纹影图所示，这一结果与文献中描述相似[10]。观察图 8空间分布，w1、w2分别为负正压力脉动，随着激波的抖振，压力脉动也随之呈现周期性运动，如图 8纹影图所示，t0至t0+2/6T时刻，负的压力脉动w1向上游移动，脉动幅值及区间增大，到t0+3/6T时刻，w1运动至翼型中部，脉动幅值及区间有所缩减，同时压力脉动也将在翼型下面进行同样的周期性运动。
|
| 图 8 翼型对称面数值纹影图与DMD第3阶模态不同时刻空间分布 Fig. 8 Numerical schlieren of symmetry plane of airfoil and spatial distribution of third-order mode of DMD at different moments |
DMD模态包含单一的流场特征频率，少数的几阶模态就能刻画流场的重要信息。图 9给出了损失函数随重构模态数目的变化曲线，Eloss为损失函数值，n为重构模态数目，损失函数定义由文献[23]给出，可以看出使用前4阶模态的7个模态(1个静模态加3对共轭模态)重构就能使损失函数降到4%以内，继续增加重构的模态数目损失函数变化较小。采用7个模态重构原始流场，与数值计算结果进行对比，流场重构的均方根误差(相对于CFD计算结果)分布见图 10，RMSE表示均方根误差值。在激波运动的区域以及尾迹区，重构的流场误差较大，取图 10中A、B、C三处作为测点，考查重构流场对压力场的捕捉能力，压力随时间的脉动如图 11所示。A点靠近翼型前缘，重构误差较小，重构的压力值与计算值基本吻合；B、C两点分别靠近翼型中部和尾缘，处于激波运动的区域内，重构流场整体上刻画了压力随时间的脉动，但在压力脉动的峰值上重构误差较大。对压力场的重构，表明DMD对于激波间断处的捕捉误差较大。
|
| 图 9 损失函数随模态数目的变化 Fig. 9 Variation of loss function with number of modes |
|
| 图 10 流场重构的均方根误差 Fig. 10 Root mean square errors of flow reconstruction |
|
| 图 11 测点压力随时间的变化 Fig. 11 Variation of observation point pressure with time |
1) DMD能准确捕捉包含流场特征频率的各阶主模态。第1阶模态的主频与激波自激振荡的频率相同，在激波的自激振荡过程占主导作用，且其正负压力脉动随激波的抖振在翼型上下面周期性运动。
2) DMD提取的前5阶模态的增长率/衰减率在0附近，对流场非定常流的发展起着持续作用，且其模态系数随时间变化呈现简谐振动特征。
3) DMD对于激波间断处的捕捉能力较弱。采用前7阶模态重构流场，能使损失函数降到4%以内，在激波运动区域重构误差较大。
| [1] |
张伟伟, 高传强, 叶正寅. 机翼跨声速抖振研究进展[J]. 航空学报, 2015, 36(4): 1056-1075. ZHANG W W, GAO C Q, YE Z Y. Research advances of wing/airfoil transic buffet[J]. Acta Aeronautica et Astronautica Sinica, 2015, 36(4): 1056-1075. (in Chinese) |
| [2] |
MCDEVITT J B, LEVY J R, DEIWERT G S. Transonic flow about a thick circular-arc airfoil[J]. AIAA Journal, 1976, 14(5): 606-613. DOI:10.2514/3.61402 |
| [3] |
TIJDEMAN H, SEEBASS R. Transonic flow past oscillating airfoils[J]. Annual Review of Fluid Mechanics, 1980, 12: 181-222. DOI:10.1146/annurev.fl.12.010180.001145 |
| [4] |
LEE B. Self-sustained shock oscillations on airfoils at transonic speeds[J]. Progress in Aerospace Sciences, 2001, 37(2): 147-196. DOI:10.1016/S0376-0421(01)00003-3 |
| [5] |
JACQUIN L, MOLTON P, DECK S, et al. Experimental study of shock oscillation over a transonic supercritical profile[J]. AIAA Journal, 2009, 47(9): 1985-1994. DOI:10.2514/1.30190 |
| [6] |
HARTMANN A, KLAAS M. Time-resolved stereo PIV measurements of shock-boundary layer interaction on a supercritical airfoil[J]. Experiments in Fluids, 2012, 52(3): 591-604. DOI:10.1007/s00348-011-1074-6 |
| [7] |
CHUNG I, LEE D, REU T.Prediction of transonic buffet onset for an airfoil with shock induced separation bubble using steady Navier-Stokes solver: AIAA-2002-2934[R].Reston: AIAA, 2002.
|
| [8] |
XIAO Q, TSAI H M, LIU F. Numerical study of transonic buffet on a supercritical airfoil[J]. AIAA Journal, 2006, 44(3): 620-628. DOI:10.2514/1.16658 |
| [9] |
XIONG J T, LIU F, LUO S J.Computation of NACA0012 airfoil transonic buffet phenomenon with unsteady Navier-Stokes equations: AIAA-2012-0699[R].Reston: AIAA, 2012.
|
| [10] |
CHEN L W, XU C Y, LU X Y. Numerical investigation of the compressible flow past an aerofoil[J]. Journal of Fluid Mechanics, 2010, 643(3): 97-126. |
| [11] |
ROWLEY C W, COLONIUS T, MURRAY R M, et al.Proper orthogonal decomposition of 2D compressible DNS of the flow over a rectangular cavity[C]//Division of Fluid Dynamics Meeting, 1999. https://www.researchgate.net/publication/241459502_Proper_Orthogonal_Decomposition_of_2D_Compressible_DNS_of_the_Flow_over_a_Rectangular_Cavity
|
| [12] |
SCHMID P J. Dynamic mode decomposition of numerical and experimental data[J]. Journal of Fluid Mechanics, 2010, 656(10): 5-28. |
| [13] |
潘翀, 陈皇, 王晋军.复杂流场的动力学模态分解[C]//第八届全国实验流体力学学术会议论文集.广州: 中国科学院南海海洋研究所, 2010: 77-82. PAN C, CHEN H, WANG J J.Dynamical mode decomposition of complex flow field[C]//8th National Conference on Experimental Fluid Mechanics.Guangzhou: South China Sea Institute of Oceanology, 2010: 77-82(in Chinese). http://www.wanfangdata.com.cn/details/detail.do?_type=conference&id=7385271 |
| [14] |
LIU H K, YAN C, ZHAO Y T, et al. Analysis of pressure fluctuation in transonic cavity flows using modal decomposition[J]. Aerospace Science & Technology, 2018, 77: 819-835. |
| [15] |
寇家庆, 张伟伟, 高传强. 基于POD和DMD方法的跨声速抖振模态分析[J]. 航空学报, 2016, 37(9): 2679-2689. KOU J Q, ZHANG W W, GAO C Q. Modal analysis of transonic buffet based on POD and DMD method[J]. Acta Aeronautica et Astronautica Sinica, 2016, 37(9): 2679-2689. (in Chinese) |
| [16] |
SPEZIALE C G, ABID R, ANDERSON E C. Critical evaluation of two-equation models for near-wall turbulence[J]. AIAA Jornal, 1992, 30(2): 324-331. DOI:10.2514/3.10922 |
| [17] |
SPALART P R, DECK S, SHUR M L, et al. A new version of detached-eddy simulation, resistant to ambiguous grid densties[J]. Theoretical and Computational Fluid Dynamics, 2006, 20: 181-195. DOI:10.1007/s00162-006-0015-0 |
| [18] |
VAN LEER B. Towards the ultimate conservative difference scheme.V.A second-order sequel to Godunov's method[J]. Journal of Computational Physics, 1979, 32(1): 101-136. |
| [19] |
YOON S, JAMESON A. Lower-upper symmetric-Gauss-Seidel method for the Euler and Navier-Stokes equations[J]. AIAA Journal, 1988, 26(9): 1025-1026. DOI:10.2514/3.10007 |
| [20] |
JAMESON A.Time dependent calculations using multigrid with applications to unsteady flows past airfoils and wings: AIAA 1991-1596[R].Reston: AIAA, 1991.
|
| [21] |
ROWLEY C W, MEZI C, BAGHERI S, et al. Spectral analysis of nonlinear flows[J]. Journal of Fluid Mechanics, 2009, 641(1): 115-127. |
| [22] |
CHEN K K, TU J H, ROWLEY C W. Variants of dynamic mode decomposition:Boundary condition, Koopman, and Fourier analyses[J]. Journal of Nonlinear Science, 2012, 22(6): 887-915. DOI:10.1007/s00332-012-9130-9 |
| [23] |
JOVANOVIC M R, SCHMID P J, NICHOLS J W. Sparsity-promoting dynamic mode decomposition[J]. Physics of Fluids, 2014, 26(2): 561-571. |