近失速形态下冰脊分离非定常流的IDDES和模态分析

doi:10.16183/j.cnki.jsjtu.2020.427

近失速形态下冰脊分离非定常流的IDDES和模态分析

谭雪^a, 张辰^,^a, 徐文浩^b, 王福新^a, 文敏华^c

a.航空航天学院,上海交通大学,上海 200240

b.船舶海洋与建筑工程学院,上海交通大学,上海 200240

c.高性能计算中心,上海交通大学,上海 200240

Unsteadiness and Modal Analysis of Ridge Ice-Induced Separation in Post-Stall Conditions via IDDES

TAN Xue^a, ZHANG Chen^,^a, XU Wenhao^b, WANG Fuxin^a, WEN Minhua^c

a. School of Aeronautics and Astronautics, Shanghai Jiao Tong University, Shanghai 200240, China

b. School of Naval Architecture, Ocean and Civil Engineering, Shanghai Jiao Tong University, Shanghai 200240, China

c. Center for High Performance Computing, Shanghai Jiao Tong University, Shanghai 200240, China

通讯作者: 张辰,男,助理研究员,电话(Tel.):021-34206388;E-mail:piressjtu@sjtu.edu.cn.

责任编辑: 石易文

收稿日期: 2020-12-21

基金资助:

结冰与防除冰重点实验室开放基金(IADL20190203)
上海交通大学新进青年教师启动计划资助项目(20X100040034)

Received: 2020-12-21

作者简介 About authors

谭雪(1996-),女,四川省德阳市人,硕士生,主要从事飞机结冰研究.

摘要

采用改进延迟脱体涡模拟 (IDDES) 方法,对近失速条件下溢流冰脊诱导的剪切层振荡现象进行高分辨率模拟,描述高雷诺数下冰脊分离流大尺度分离的流场演化特征.研究表明,近失速形态下,冰脊和下翼面尾缘同时诱导出剪切流动,冰脊诱导的剪切层并未再附到壁面,与下翼面上洗流动相互干扰,形成大尺度低能态结构.结合频谱分析进一步发现,剪切层内的压力脉动存在两种典型的脉动频率,与Kelvin-Helmholtz不稳定性相关,具体表现为涡配对和涡脱落.基于正交分解得到的压力脉动主导模态为剪切层之间的大尺度相干结构.主导模态的时间系数与升力系数的功率密度谱峰值保持一致,这说明尾缘区大尺度相干结构与升力波动具有相关性.

关键词： 溢流冰脊; 剪切层振荡; 数值模拟; 频谱分析; 模态分析

Abstract

High-resolution simulation of shear layer oscillation induced by ridge ice in post-stall condition is conducted via the improved delayed detached-eddy simulation (IDDES) method. The flow-field evolution characteristics of large scale separation in high Reynolds number condition are described. It is shown that the ridge ice and trailing edge of the lower surface induce the development of shear flow at the same time. The wall is not reattached by the shear layer induced by ridge ice, and the “up-wash” flow from the lower surface is interacted with the shear layer, which lead to the formation of large-scale coherent structures. Combined with the spectral analysis, the pressure pulsation located in the shear layer is characterized by two typical frequencies, which are associated with Kelvin-Helmholtz instability and appear as the vortex pairing and shedding. Based on the proper orthogonal decomposition, the dominant mode of pressure pulsation between shear layers is extracted as large-scale coherent structures. The same peak value is shown in power density spectrum of dominant mode temporal coefficient and lift coefficient, which indicates that the large-scale coherent structure is connected with lift fluctuation.

Keywords： ridge ice; shear layer oscillation; numerical simulation; spectrum analysis; modal analysis

PDF (15798KB) 元数据多维度评价相关文章导出 EndNote| Ris| Bibtex 收藏本文

本文引用格式

谭雪, 张辰, 徐文浩, 王福新, 文敏华. 近失速形态下冰脊分离非定常流的IDDES和模态分析[J]. 上海交通大学学报（自然版）, 2021, 55(11): 1333-1342 doi:10.16183/j.cnki.jsjtu.2020.427

TAN Xue, ZHANG Chen, XU Wenhao, WANG Fuxin, WEN Minhua. Unsteadiness and Modal Analysis of Ridge Ice-Induced Separation in Post-Stall Conditions via IDDES[J]. Journal of shanghai Jiaotong University, 2021, 55(11): 1333-1342 doi:10.16183/j.cnki.jsjtu.2020.427

溢流冰脊通常形成于机翼防除冰装置后方^[1],其引起的流动大尺度分离容易导致飞机控制面失效,是过冷大水滴环境影响飞机飞行安全的典型特征之一^[2].冰脊虽然结冰量较少,对气动性能产生的影响却十分严重^[3]——冰脊前发生流动分离,并诱导剪切层振荡,伴随着复杂的涡脱落、破碎和再附着现象,是研究人员关注的热点问题^[4,5].美国联邦航空管理局(FAA)在最新发布的结冰适航修正案^[6]和咨询通告中多次强调:溢流冰比传统认知的临界冰型更加危险,其造成的压力低频脉动可能是导致某些控制面失效,并最终引发飞行事故的根本原因.

在早期研究中,冰脊所造成的气动损失被重点关注,对分离流的探讨主要集中于最大升力系数和失速攻角的损失,其中的代表性成果是Bragg团队基于伊利诺伊大学低湍流度风洞试验的相关工作.文献[7]发现NACA23012翼型的最大升力系数从1.47下降至0.7附近,下降了52%.文献[8]结合压力系数分布,进一步指出雷诺数对冰脊翼型近失速状态下的升阻力影响并不大,这主要归因于冰脊翼型的前缘发生大尺度分离,减少了壁面流动对雷诺数的依赖.文献[9]后续进一步研究了带舵面翼型发生溢流冰现象的气动损失,证明了这类冰型的破坏性.

为了解析这一物理现象,人们尝试引入计算流体动力学(CFD)方法对非定常流动过程作进一步研究.由于雷诺平均Navier-Stoks(RANS)方程的时均处理方式难以描述分离流的不稳定性,文献[10]把RANS方程和大涡模拟(LES)的混合算法引入冰脊翼型的非定常流计算,最早实现了对冰脊剪切层发展、振荡、脱落和附着过程的模拟.研究表明,冰脊翼型存在瞬态气动力大幅波动的现象.文献[11]提出了改进延迟脱体涡模拟(IDDES)的模拟方案,用于改进翼面压力分布的模拟.文献[12]提出采用壁面模化大涡模拟(WMLES)改善网格计算不足时的模化应力导致的计算“灰区”问题.文献[13,14]分别提出5阶加权基本无振荡 (WENO) 低耗散格式、自适应长度尺度等方法进一步改善剪切层失稳破碎的模拟能力.这些研究均为冰脊剪切流条件下的高分辨率解析提供了重要的技术支撑.

在最新的试验研究中,文献[15]发现冰脊流场中存在一种区别于常规模态的低频模态,表现为剪切层挥舞形态,与压力波动存在关联.剪切区的谱分析结果表明,这种振荡特性的Strouhal数在0.02附近^[16],与低雷诺数下后台阶流的流动特性比较相似.结合后台阶流理论^[17]不难发现,剪切层尾迹中的涡破裂,引发流体向外卷吸速度和向内的回注速度同时发生变化,这种“流量夹带”和“卷吸流量回注”之间的瞬时不平衡造成了剪切层发展的不稳定,剪切层在这种作用下会伴随尾流发生低频自维持振荡,然而这种现象的本质并没有完全得到解释.

本文的研究目的在于采用IDDES方法^[18],对溢流冰脊诱导剪切层振荡的非定常流特性进行高分辨率模拟,期望通过结合正交分解方法,进一步获取近失速形态下,造成冰脊分离流压力波动的主导模态,深入探讨冰脊分离流动不稳定性的时域和频域特征,探索挖掘大尺度相干结构与升力波动的关联性.本文的研究揭示了近失速形态下冰脊分离流导致升力剧烈振荡的物理机制.

1 研究方法

1.1 计算条件

采用数值模拟方法对近失速形态下的冰脊翼型分离流模态特征进行计算分析.来流马赫数Ma设置为0.21,来流迎角α为5°,弦长参考雷诺数Re_c为2.0×10⁶.翼型模型采用NACA23012,冰脊在x方向位置与翼型弦长c之比x/c为0.1,冰脊高度k与翼型弦长c之比为 0.0139.该模型几何外形规则无毛刺,便于生成高质量的贴体网络,通过合适的网格分布设计,可以保证冰脊附近的网格正交性.由于文献[7]在早期的NACA23012缩比模型实验中已经证明:冰脊采用后向1/4圆时会产生最不利气动危害,并以此开展了详细的流场特征研究,得到一系列气动数据,所以本研究选择该冰脊外形特征进行对比.计算条件与文献[10]分离涡模拟(DES)的模拟结果 $(Ma = 0.21, Re = 2.0 \times 10^{6})$ 相同.

1.2 计算方法

计算模型的几何外形如图1所示,其中: $y / c 为 y$ 方向位置与弦长之比.采用基于CFL3D框架开发的高分辨率流场数值求解器^[19]实现三维冰脊流场的计算模拟.基于连续性方程、动量方程和理想气体方程实现ρ、p、u、v、w、T这6个未知量的计算封闭,其中:ρ为密度;p为压力;u,v,w分别为x,y,z方向的速度;T为温度.忽略体积力和外部热源,其守恒形式可写成矢量形式如下:

(1)

\frac{\partial H}{\partial t} + \frac{\partial F_{i}}{\partial x_{i}} = \frac{\partial F_{v, i}}{\partial x_{i}}

(2)

\begin{array}{l} \begin{matrix} \begin{array}{l} H = {[\begin{array}{l} ρ & ρ u_{1} & ρ u_{2} & ρ u_{3} & E \end{array}]}^{T} \\ F_{i} = [\begin{array}{l} ρ u_{i} \\ ρ u_{i} u_{1} + p δ_{i 1} \\ ρ u_{i} u_{2} + p δ_{i 2} \\ ρ u_{i} u_{3} + p δ_{i 3} \\ (E + p) u_{i} \end{array}] \\ F_{v, i} = {[\begin{array}{l} 0 & τ_{i 1} & τ_{i 2} & τ_{i 3} & W_{i} \end{array}]}^{T} \end{array}\} \end{matrix} \\ i = 1,2, 3 \end{array}

式中: H为状态变量;F_i 为i方向的对流矢量通量;F_v,i 为i方向的黏性矢量通量;i=1,2,3分别为笛卡尔坐标系下的x,y,z方向;x_i 为笛卡尔坐标系下的不同方向;t为时间;E为能量;u_i 为笛卡尔坐标系下x,y,z方向的速度分量;δ_ij 为克罗内克尔符号,i≠j,j=1,2,3;W_i 为i方向黏性应力所做的功.按照Stokes假设,黏性应力τ_ij 及其所做的功W_i可以表示为

(3)

τ_{ij} = μ (\frac{\partial u_{i}}{\partial x_{j}} + \frac{\partial u_{j}}{\partial x_{i}}) - \frac{2}{3} μ \frac{\partial u_{k}}{\partial x_{k}} δ_{ij}

(4)

W_{i} = u_{j} τ_{ij} + κ \frac{\partial T}{\partial x_{i}}

式中: x_j 为笛卡尔坐标系下的不同方向,i≠j,j=1,2,3;u_k 为笛卡尔坐标系下x,y,z方向的速度分量且k=1,2,3,i≠j≠k;u_j 为在笛卡尔坐标系下x,y,z方向的速度分量,i≠j;μ为动力黏度;κ为导热系数.对于理想气体,根据状态方程,对单位体积气体总能量定义的数学表达式做恒等变换后可以得到:

(5)

\begin{matrix} p = (γ - 1) [E - \frac{1}{2} ρ (u_{1}^{2} + u_{2}^{2} + u_{3}^{2})] \end{matrix}

式中: γ为比热比.动力黏度μ与温度T有关,可通过下式计算,则有:

(6)

\begin{matrix} \frac{μ}{μ_{0}} \approx {(\frac{T}{T_{0}})}^{1.5} (\frac{T_{0} + T_{s}}{T + T_{s}}) \end{matrix}

式中:T₀=273.16 K,为参考温度;对于空气,有μ₀=17.161 μPa·s,为参考温度对应的动力黏度;T_s=124 K,为与空气对应的常数.

图1

新窗口打开| 下载原图ZIP| 生成PPT

图1 NACA23012-QR溢流冰脊翼型的几何形状

Fig.1 Geometry of NACA23012-QR ridge ice airfoil

湍流涡黏性采用重新定义的长度尺度和LES/RANS特征自适应混合函数对Spalart-Allmaras(S-A)模型进行模化.LES/RANS特征自适应混合函数l_hyb可由下式定义:

(7)

\begin{matrix} l_{hyb} = {\tilde{f}}_{d} (1 + f_{e}) l_{RANS} + (1 - {\tilde{f}}_{d}) l_{LES} \end{matrix}

式中:l_RANS为RANS的长度尺度;l_LES为LES的长度尺度;f_e为升降函数; ${\tilde{f}}_{d}$ 为混合函数,

(8)

\begin{matrix} {\tilde{f}}_{d} = \max \{(1 - f_{dt}), f_{B}\} \end{matrix}

f_dt为与湍流模型相关的函数,f_B为经验混合函数.当进口条件包含湍流时, ${\tilde{f}}_{d}$ =f_B.此时,亚格子长度尺度变为WMLES的长度尺度,即l_hyb=l_WM,则有:

(9)

\begin{matrix} l_{WM} = f_{B} (1 + f_{e}) l_{RANS} + (1 - f_{B}) l_{LES} \end{matrix}

此时,WMLES被激活,并快速从RANS模式切换到LES模式,确保分离区除近壁面区域外的大部分湍流流动能够计算准确.

当进口条件不包含湍流时,亚格子长度尺度变为DDES的长度尺度,即l_hyb=l_DDES,

(10)

\begin{matrix} l_{DDES} = {\tilde{f}}_{d} l_{RANS} + (1 - {\tilde{f}}_{d}) l_{LES} \end{matrix}

这时,亚格子长度尺度与延迟脱体涡(DDES)方法保持一致,保留了原始DDES方法模拟分离流场的能力.

综上所述,新的亚格子尺度可以用下式表达:

(11)

Δ = \min \{\max \{C_{w} d_{w}, C_{w} h_{\max}, h_{wn}\}, h_{\max}\}

式中:h_wn为沿壁面的法向单元网格尺度;C_w为壁面常数;h_max为网格在各方向上的最大尺度;d_w为网格距离壁面的长度.由式(11)可以看出,新定义的长度尺度中既考虑了局部网格的影响,又考虑了壁面距离的影响.近壁面的亚格子长度尺度减小,导致近壁面涡黏系数急剧下降.这种效果增强了长度尺度对近壁面和自由湍流模拟的适应性,更加适合用于模拟复杂的湍流分离流动.

数值离散方面,空间上对流项计算采用基于5阶WENO和6阶中心差分的混合格式进行离散^[19],黏性项采用4阶中心格式;时间上采用双时间步进策略^[20]和隐式近似因子分解内迭代法以2阶精度推进求解.物理时间导数项采用2阶精度的三点向后差分格式离散,伪时间导数项采用1阶向后差分格式逼近.

1.3 计算网格

空间离散采用结构化网格.展向计算域长度为10%翼型弦长,展向面采用周期边界条件进行关联.网格划分的分区策略如图2所示.由图2可知,翼型上方冰脊附近和尾迹区域被重点加密,以保证对冰脊诱导的非定常分离流动具备足够的分辨率;这两个区域的网格点被精心布置,以保证其具有绝对正交性;且通过局部空间的网格加密可以有效改善DES方法受灰区的影响范围;近壁区域网格按照无量纲壁面距离y⁺准则生成(第1层网格到壁面的距离约为3×10^-6),以保证壁面压力分布的预测要求.

图2

新窗口打开| 下载原图ZIP| 生成PPT

图2 冰脊翼型网格策略示意图

Fig.2 Schematic diagram of mesh strategy for ridge ice airfoil

为了验证网格无关性,分别取粗网格、中网格以及细网格进行对比计算.3套计算网格的详细信息如表1所示, 其中:Δy₁为第1层网格质心到近壁面的距离; N为翼面上布置的网格点数量;ξ为增长率;Δt^*为无量纲时间步长; ${\bar{C}}_{L}$ 为平均升力系数.3套网格的总网格量约为8×10⁶、1.1×10⁷以及 1.3×10⁷.由表1可知,中网格和密网格的平均升力系数结果已基本趋于一致,且与实验结果基本符合, 因此这里选用中网格开展后续研究.

表1 无关性验证的网格参数

Tab.1 Irrelevance verification of mesh parameter

网格	网格量	Δy₁	N	ξ	Δt^*	${\bar{C}}_{L}$
粗网格	8.0×10⁶	3×10^-6	535	1.1	0.01/340	0.369
中网格	1.1×10⁷	3×10^-6	571	1.05	0.01/340	0.3015
细网格	1.3×10⁷	3×10^-6	602	1.05	0.01/340	0.3011
实验	—	—	—	—	—	0.3010

新窗口打开| 下载CSV

1.4 计算验证

在α=5° 情况下,NACA23012-QR翼型的压力系数分布如图3(a)所示,其中:C_p为上翼面的压力系数.由图3(a)可以看出,C_p计算值在x/c=0~0.45位置与实验结果的符合性较好;在x/c=0.45~1位置比试验值略微偏小. IDDES获得的下翼面压力分布要明显优于DES,整个近失速形态下的压力分布与试验测量的相差不大.冰脊翼型的计算平均升力线 ${\bar{C}}_{L}$ -α特征如图3(b)所示.从图3(b)中可以看出,IDDES方法准确地预测了翼型的平均升力系数,计算结果有效.

图3

新窗口打开| 下载原图ZIP| 生成PPT

图3 冰脊翼型的C_p及 ${\bar{C}}_{L}$ -α验证

Fig.3 Validation with C_p and ${\bar{C}}_{L}$ -α of ridge ice airfoil

2 结果与讨论

2.1 流场分析

基于Q准则的NACA23012-QR冰脊翼型瞬时流场特征(Re_c=2.0×10⁶,Ma=0.21,α=5° ),如图4所示. 其中:Q=0.01,为流场中速度梯度张量的第2矩阵不变;涡结构采用马赫数进行着色.从图4中可以清晰地看出冰脊分离流的空间演化特征,相比于传统的角状冰^[12],冰脊在前缘提前诱导出流动分离,在冰脊前方产生驻涡并诱导出大速度的剪切层流动.这种流动呈现出剪切层Kelvin-Helmholtz不稳定性演化出的管状发卡涡结构,向空间各向扭曲,最终演化为湍流漩涡,与角状冰的涡脱落和破碎现象有一定相似之处.由于前方驻涡的影响,冰脊流的滞止线与翼型切面的夹角增大,使得分离流相比角冰更难以完全再附到翼型表面,这很可能是其在较小攻角下发生尾流剧烈不稳定振荡的主要原因.

图4

新窗口打开| 下载原图ZIP| 生成PPT

图4 基于马赫数着色的溢流冰脊Q准则瞬时流场(0≤Ma≤0.32)

Fig.4 Instantaneous flow of ridge ice Q criterion colored by Mach number (0≤Ma≤0.32)

在α=5° 情况下,冰脊分离流平均马赫数云图如图5所示.由图5可以看出,冰脊和下翼面尾缘同时诱导出两个沿流向的剪切层.其中,冰脊诱导的剪切层并未再附到壁面,与下翼面的剪切流动在远场汇聚.两个剪切流之间形成大尺度的低能态结构,导致整个上翼面几乎不存在附面层流动.

图5

新窗口打开| 下载原图ZIP| 生成PPT

图5 时均马赫数云图(0≤Ma≤0.28)

Fig.5 Contours of time-averaged Mach number (0≤Ma≤0.28)

为了研究其不稳定性,进一步分析其湍动能E_k.如图6所示.由图6可知,在靠近冰脊的回流区的湍动能较弱,这说明剪切层在x/c=0.1~0.3区域的流量夹带很少影响回流区.在x/c=0.4以后,剪切层充分发展,翼型上方的湍动能明显增强.湍动能的最大值出现在翼型上方的后缘处,这可能是由于冰脊和下翼面尾缘的剪切层发生相互干扰,从而对整个流场的不稳定性产生影响.关于这种影响的量化描述,在下文中结合主导模态做更详细的分析.

图6

新窗口打开| 下载原图ZIP| 生成PPT

图6 湍动能云图(0≤E_k≤0.005,从暗到亮)

Fig.6 Contours of turbulent kinetic energy (0≤E_k≤0.005,from darkness to brightness)

升力系数的历史曲线图如图7(a)所示,其中:Φ为升力系数功率密度谱(PSD).由图7(a)可知,升力波动的范围在0.2~0.6之间,与Loth采用DES方法计算的升力系数范围大致吻合.利用Fourier分析获得的Φ如图7(b)所示, 其中:频率用斯坦顿数St进行归一化.由图7(b)可知,冰脊翼型的瞬态升力系数虽然出现了多个频率,但主导频率出现在0.02~0.025之间.

图7

新窗口打开| 下载原图ZIP| 生成PPT

图7 升力系数历史曲线及其PSD

Fig.7 History of lift coefficient and its PSD

压力脉动p'的云图如图8所示.由图8可知,压力脉动在x/c=0.4达到最大,这可能与此处剪切层涡对流形正在逐渐演化为大尺度涡结构相关.在x/c=0.2,0.3,0.4,0.5,0.8,1.0处(L1~L6)提取的压力系数脉动C_p,rms特征分布如图9所示,其中:rms为均方根值(见图8中L1~L6). 由图9可知, x/c=0.2~0.3处的压力脉动受剪切层约束,压力脉动空间主要位于近壁区;x/c=0.4后方,由于在尾迹中的大尺度结构相互干扰,压力脉动分布趋于平缓.

图8

新窗口打开| 下载原图ZIP| 生成PPT

图8 压力脉动云图

Fig.8 Contour of pressure pulsation

图9

新窗口打开| 下载原图ZIP| 生成PPT

图9 不同位置的C_p,rms分布

Fig.9 Distribution of C_p,rms at different locations

2.2 频谱分析

t₀,t₁,t₂ 时刻的瞬时压力脉动云图如图10所示.每个时刻间隔1个时间步长,从t₀~t₂分离涡逐渐向下游移动,前序分离涡在发展中逐渐变大.图11提取了P1,P2,P3,P4这4个位置的压力脉动能谱.P1点的压力脉动能谱存在St=0.1325及St=0.025 两个峰值;文献[12]在研究角冰时将其归因于剪切层的涡配对和脱落,这在本研究中同样成立.P2和P3点位置只存在St=0.025的峰值,说明此时剪切层已演化为涡脱落;P4点位于剪切层尾迹中,只能检测到宽频信号,此时剪切层已经演化为低能态大尺度结构.

图10

新窗口打开| 下载原图ZIP| 生成PPT

图10 压力脉动的时序演化

Fig.10 Temporal evolution of pressure pulsation

图11

新窗口打开| 下载原图ZIP| 生成PPT

图11 P1~P4的压力脉动功率密度谱图

Fig.11 PSD of pressure fluctuations at locations of P1 to P4

2.3 模态分析

为了进一步分析冰脊剪切流的不稳定性特征,采用快照正交分解(POD)方法^[21]对流场的压力脉动时域信息进行处理.POD的基本思想是将时间和空间上连续的物理量分解.给定M个时刻的流场快照,任意脉动量可以写作M维POD基和对应模态系数的乘积,则有:

(12)

\begin{matrix} p' (x, y, t_{l}) = \overset{M}{\sum_{n = 1}} a_{n} (t_{l}) φ_{n} (x, y) \end{matrix}

式中: p'(x,y,t_l) 为第t_l个时刻流场的压力脉动值;φ_n (x,y)为第n 个POD模态基函数; a_n (t_l)为第t_l 个时刻中第n个POD基的模态系数; l,n=1,2,…,M,分别为原始流场的第l个时刻以及POD处理后的第n个POD模态.要计算POD基,首先需要定义相关性矩阵R,具体表示为

(13)

\begin{matrix} R = V^{T} V \end{matrix}

式中:

V=[p'(x, y, t₁) p'(x, y, t₂) … p'(x, y, t_M)],

为减掉平均流场的脉动量组成的序列向量.

由于R为Gram矩阵,其特征值必大于等于0,特征值对应的特征向量正交.在POD方法中,相关性矩阵R的特征向量决定了POD模态的形式,而R的特征值大小表征对应模态的能量.对于矩阵R_N_×_N, 可计算获得M个特征值和特征向量,即得到M个POD模态.由于特征值越大,对应于权重更高的POD模态,只需保留少数高能量的模态,就能有效捕捉原始数据中的主要特征.对相关性矩阵作奇异值分解,可以得到:

(14)

\begin{matrix} R A^{n} = λ_{n} A^{n} \end{matrix}

式中:Aⁿ=[a_n(t₁) a_n(t₂) … a_n(t_M)]^T,为模态系数矩阵;λ_n为特征值.POD的基表达式如下:

(15)

\begin{matrix} φ_{n} (x, y) = \frac{1}{M λ_{n}} \overset{M}{\sum_{l = 1}} A_{l}^{n} p' (x, y, t_{l}) \end{matrix}

由式(6)~(9)可知,POD方法能够将流场脉动分解为3个特征量:① POD特征值;② POD基模态;③ 各模态对应的时间系数.

基于奈奎斯特取样定理的样本数验证如图12(a)所示.通过对比200,400,800张样本量的模态能量占比 $λ_{n} / \sum λ_{n},$ 基本排除了快照样本数对分析结果的影响,后续研究采用800×800的取样标准.图12(b)为POD前20阶模态能量占比以及前100阶模态累计能量 $\sum λ_{j}$ 如图12(b)所示.由图12(b)可以看出,低阶模态的特征值较大,是流场不稳定性的主要因素.前20阶模态包含了流场中80%的压力脉动能力,前200阶模态占据了99%,因此本研究主要关注其前200阶的POD模态特征.

图12

新窗口打开| 下载原图ZIP| 生成PPT

图12 POD模态归一化能量及其累计能量

Fig.12 Normalized energy and its cumulative distribution of POD mode

在α=5° 情况下,流场脉动的POD基模态如图13所示.模态1和模态2的能量占比分别为18.2%和13.3%,其空间分布呈现正负间隔的能量脉动序列,主要分布在翼型后缘处, 表现为明显的大尺度相干结构,影响区域涉及后缘的近壁区,对后缘壁面压力系数造成影响.模态3和模态4的空间分布与前两阶模态特点相似,尺度略小但更靠近壁面.模态49和模态50的能量占比为0.21%,且其影响区间范围与近壁面的相关性较小.模态199和200的能量占比仅为0.09%,其压力脉动范围几乎不涉及壁面,可推断其对壁面压力波动几乎不产生影响.

图13

新窗口打开| 下载原图ZIP| 生成PPT

图13 冰脊分离流的典型POD基模态

Fig.13 Typical POD modes of ice-ridge separation flow

图14对比了升力系数频谱与1/2阶POD模态、3/4阶POD模态、49/50阶POD模态以及199/200阶POD模态的时间系数功率密度谱,其中:a_n(t)为第n个POD基的模态时间系数.由图14可知,升力脉动的频谱特征几乎与前4阶模态的POD时间系数重合,这充分说明大尺度相干结构对升力脉动的影响.相比较而言,高阶模态的时间系数功率谱与升力系数频谱并不吻合,这暗示了在近失速形态下,冰脊诱导的剪切层振荡结构对升力波动的影响并不大.

图14

新窗口打开| 下载原图ZIP| 生成PPT

图14 典型模态时间系数与升力系数能谱图的比较

Fig.14 Comparison of typical mode temporal coefficients and lift coefficients

3 结论

本文采用改进延迟脱体涡模拟方法,对近失速形态下,溢流冰脊诱导剪切层振荡的非定常流特性进行高分辨率模拟.结合湍流统计、谱分析和正交分解等方法,本研究进一步分析了冰脊分离流的非定常特性,提取了剪切区的压力振荡频谱特征和全域流场压力脉动的主导模态,探索挖掘了POD基模态与升力波动的关联性.本研究主要获得以下结论.

(1) 近失速形态下,溢流冰脊和下翼面尾缘同时诱导出两个剪切流动.冰脊诱导的剪切层与下翼面上洗流动相互干扰,形成大尺度低能态结构,从而导致整个上翼面几乎不存在附面层流动.

(2) 冰脊诱导的剪切层中存在着两种典型的脉动频率与Kelvin-Helmholtz不稳定性相关,具体表现为剪切层诱导的涡配对和涡脱落;剪切层在演化为大尺度结构位置x/c=0.4的压力脉动最为剧烈;x/c=0.4之前受剪切层约束局限于近壁区,x/c=0.4之后受大尺度结构干扰趋于平缓.

(3) 基于正交分解提取获得的主导模态空间形态为冰脊剪切流尾缘区相互作用的大尺度结构,其时间系数谱与升力脉动的功率密度谱的峰值高度相似,这表明真正造成冰脊翼型升力波动的主要因素来自于尾缘中相互作用的大尺度相干结构.

参考文献

原文顺序

文献年度倒序

文中引用次数倒序

被引期刊影响因子

[1]

顾洪宇, 桑为民, 庞润, 等.

机翼热气防冰及冰脊形成数值模拟

[J]. 气体物理, 2019, 4(4):41-49.

[本文引用: 1]

Hongyu

, SANG

Weimin

, PANG

Run

, et al.

Numerical simulation of wing hot air anti-icing and ice ridge formation

[J]. Physics of Gases, 2019, 4(4):41-49.

[本文引用: 1]

[2]

ZHANG

, LIU

Effect of drop size on the impact thermodynamics for supercooled large droplet in aircraft icing

[J]. Physics of Fluids, 2016, 28(6):062107.

DOI:10.1063/1.4953411 URL [本文引用: 1]

[3]

POURYOUSSEFI S

, MIRZAEI

, NAZEMI M

, et al.

Experimental study of ice accretion effects on aerodynamic performance of an NACA 23012 airfoil

[J]. Chinese Journal of Aeronautics, 2016, 29(3):585-595.

DOI:10.1016/j.cja.2016.03.002 URL [本文引用: 1]

[4]

ZHANG

, HABASHI W

, KHURRAM R

Zonal detached-eddy simulation of turbulent unsteady flow over iced airfoils

[J]. Journal of Aircraft, 2015, 53(1):168-181.

DOI:10.2514/1.C033253 URL [本文引用: 1]

[5]

李焱鑫, 顾新.

民机SLD结冰研究和适航验证的发展与挑战

[J]. 中国民航大学学报, 2020, 38(4):48-53.

[本文引用: 1]

Yanxin

, GU

Xin

Development and challenges of supercooled large drop icing research for civil airplane airworthiness certification

[J]. Journal of Civil Aviation University of China, 2020, 38(4):48-53.

[本文引用: 1]

[6]

Federal Aviation Administration. Federal Aviation Administration on airplane and engine certification requirements in supercooled large drop, mixed phase, and ice crystal icing condition [EB/OL]. (2014-04-11)[2019-12-13]. http://www.federalregister.gov/documents/2014/11/04/2014-25789/airplane-and-engine-certification-requirements-in-supercooled-large drop-mixed-phase-and-ice-crystal-icing-condition.

URL [本文引用: 1]

[7]

BROEREN

, BRAGG

Effect of airfoil geometry on performance with simulated intercycle ice accretions

[C]// 41st Aerospace Sciences Meeting and Exhibit. Reston, Virginia, USA: AIAA, 2003: 1-16.

[本文引用: 2]

[8]

LEE

, BRAGG M

Investigation of factors affecting iced-airfoil aerodynamics

[J]. Journal of Aircraft, 2003, 40(3):499-508.

DOI:10.2514/2.3123 URL [本文引用: 1]

[9]

李冬, 张辰, 王福新, 等.

结冰对带舵面翼型流场的影响及其气动参数分析

[J]. 上海交通大学学报, 2017, 51(3):367-373.

[本文引用: 1]

Dong

, ZHANG

Chen

, WANG

Fuxin

, et al.

Effect of icing on airfoil with control surface and analysis of aerodynamic parameters

[J]. Journal of Shanghai Jiao Tong University, 2017, 51(3):367-373.

[本文引用: 1]

[10]

PAN J

, LOTH

Detached eddy simulations for iced airfoils

[J]. Journal of Aircraft, 2005, 42(6):1452-1461.

DOI:10.2514/1.11860 URL [本文引用: 2]

[11]

TRAVIN A

, SHUR M

, SPALART P

, et al. Improvement of delayed detached-eddy simulation for LES with wall modelling[EB/OL]. (2006-09-07) [2020-10-21]. https://repository.tudelft.nl/islandora/object/uuid: dbcd7d38-60cb-4107-a6cf-9f4b5dd637a9?collection=research.

URL [本文引用: 1]

[12]

XIAO M

, ZHANG Y

, ZHOU

Numerical study of iced airfoils with horn features using large-eddy simulation

[J]. Journal of Aircraft, 2018, 56(1):94-107.

DOI:10.2514/1.C034986 URL [本文引用: 3]

[13]

张恒, 李杰, 龚志斌.

基于IDDES方法的翼型结冰失速分离流动数值模拟

[J]. 空气动力学学报, 2016, 34(3):283-288.

[本文引用: 1]

ZHANG

Heng

, LI

Jie

, GONG

Zhibin

Numerical simulation of the stall separated flow around an iced airfoil based on IDDES

[J]. Acta Aerodynamica Sinica, 2016, 34(3):283-288.

[本文引用: 1]

[14]

HU S

, ZHANG

, LIU

, et al.

Study on vortex shedding mode on the wake of horn/ridge ice contamination under high-Reynolds conditions

[J]. Proceedings of the Institution of Mechanical Engineers, Part G: Journal of Aerospace Engineering, 2019, 233(13):5045-5056.

DOI:10.1177/0954410019835971 URL [本文引用: 1]

[15]

ANSELL P

, BRAGG M

Unsteady modes in flowfield about airfoil with horn-ice shape

[J]. Journal of Aircraft, 2015, 53(2):475-486.

DOI:10.2514/1.C033421 URL [本文引用: 1]

[16]

ANSELL P

, BRAGG M

Characterization of low-frequency oscillations in the flowfield about an iced airfoil

[J]. AIAA Journal, 2014, 53(3):629-637.

DOI:10.2514/1.J053206 URL [本文引用: 1]

[17]

MA X

, SCHRÖDER

Analysis of flapping motion of reattaching shear layer behind a two-dimensional backward-facing step

[J]. Physics of Fluids, 2017, 29(11):115104.

DOI:10.1063/1.4996622 URL [本文引用: 1]

[18]

HU S

, ZHANG

, LIU

, et al.

IDDES simulation of flow separation on an 3-D NACA23012 airfoil with spanwise ridge ice

[C]// 2018 Atmospheric and Space Environments Conference. Atlanta, Georgia, USA: AIAA, 2018: 1-10.

[本文引用: 1]

[19]

LI G

, FU

, WANG F

High-resolution multi-code implementation of unsteady Navier-Stokes flow solver based on paralleled overset adaptive mesh refinement and high-order low-dissipation hybrid schemes

[J]. International Journal of Computational Fluid Dynamics, 2017, 31(9):379-395.

DOI:10.1080/10618562.2017.1387251 URL [本文引用: 2]

[20]

JAMESON

Time dependent calculations using multigrid, with applications to unsteady flows past airfoils and wings

[C]// 10th Computational Fluid Dynamics Conference. Honolulu, HI, USA: AIAA, 1991: 1-14.

[本文引用: 1]

[21]

ZHANG

, HU S

, LI Y

, et al.

Modal analysis of 3-D iced-airfoil aerodynamics based on proper orthogonal decomposition techniques

[C]// AIAA Scitech 2019 Forum. San Diego, California, USA: AIAA, 2019: 1-11.

[本文引用: 1]

机翼热气防冰及冰脊形成数值模拟

2019

... 溢流冰脊通常形成于机翼防除冰装置后方^[1],其引起的流动大尺度分离容易导致飞机控制面失效,是过冷大水滴环境影响飞机飞行安全的典型特征之一^[2].冰脊虽然结冰量较少,对气动性能产生的影响却十分严重^[3]——冰脊前发生流动分离,并诱导剪切层振荡,伴随着复杂的涡脱落、破碎和再附着现象,是研究人员关注的热点问题^[4,5].美国联邦航空管理局(FAA)在最新发布的结冰适航修正案^[6]和咨询通告中多次强调:溢流冰比传统认知的临界冰型更加危险,其造成的压力低频脉动可能是导致某些控制面失效,并最终引发飞行事故的根本原因. ...

机翼热气防冰及冰脊形成数值模拟

2019

Effect of drop size on the impact thermodynamics for supercooled large droplet in aircraft icing

2016

Experimental study of ice accretion effects on aerodynamic performance of an NACA 23012 airfoil

2016

Zonal detached-eddy simulation of turbulent unsteady flow over iced airfoils

2015

民机SLD结冰研究和适航验证的发展与挑战

2020

民机SLD结冰研究和适航验证的发展与挑战

2020

Effect of airfoil geometry on performance with simulated intercycle ice accretions

2003

... 在早期研究中,冰脊所造成的气动损失被重点关注,对分离流的探讨主要集中于最大升力系数和失速攻角的损失,其中的代表性成果是Bragg团队基于伊利诺伊大学低湍流度风洞试验的相关工作.文献[7]发现NACA23012翼型的最大升力系数从1.47下降至0.7附近,下降了52%.文献[8]结合压力系数分布,进一步指出雷诺数对冰脊翼型近失速状态下的升阻力影响并不大,这主要归因于冰脊翼型的前缘发生大尺度分离,减少了壁面流动对雷诺数的依赖.文献[9]后续进一步研究了带舵面翼型发生溢流冰现象的气动损失,证明了这类冰型的破坏性. ...

... 采用数值模拟方法对近失速形态下的冰脊翼型分离流模态特征进行计算分析.来流马赫数Ma设置为0.21,来流迎角α为5°,弦长参考雷诺数Re_c为2.0×10⁶.翼型模型采用NACA23012,冰脊在x方向位置与翼型弦长c之比x/c为0.1,冰脊高度k与翼型弦长c之比为 0.0139.该模型几何外形规则无毛刺,便于生成高质量的贴体网络,通过合适的网格分布设计,可以保证冰脊附近的网格正交性.由于文献[7]在早期的NACA23012缩比模型实验中已经证明:冰脊采用后向1/4圆时会产生最不利气动危害,并以此开展了详细的流场特征研究,得到一系列气动数据,所以本研究选择该冰脊外形特征进行对比.计算条件与文献[10]分离涡模拟(DES)的模拟结果

(Ma = 0.21, Re = 2.0 \times 10^{6})

相同. ...

Investigation of factors affecting iced-airfoil aerodynamics

2003

结冰对带舵面翼型流场的影响及其气动参数分析

2017

结冰对带舵面翼型流场的影响及其气动参数分析

2017

Detached eddy simulations for iced airfoils

2005

... 为了解析这一物理现象,人们尝试引入计算流体动力学(CFD)方法对非定常流动过程作进一步研究.由于雷诺平均Navier-Stoks(RANS)方程的时均处理方式难以描述分离流的不稳定性,文献[10]把RANS方程和大涡模拟(LES)的混合算法引入冰脊翼型的非定常流计算,最早实现了对冰脊剪切层发展、振荡、脱落和附着过程的模拟.研究表明,冰脊翼型存在瞬态气动力大幅波动的现象.文献[11]提出了改进延迟脱体涡模拟(IDDES)的模拟方案,用于改进翼面压力分布的模拟.文献[12]提出采用壁面模化大涡模拟(WMLES)改善网格计算不足时的模化应力导致的计算“灰区”问题.文献[13,14]分别提出5阶加权基本无振荡 (WENO) 低耗散格式、自适应长度尺度等方法进一步改善剪切层失稳破碎的模拟能力.这些研究均为冰脊剪切流条件下的高分辨率解析提供了重要的技术支撑. ...

(Ma = 0.21, Re = 2.0 \times 10^{6})

相同. ...

Numerical study of iced airfoils with horn features using large-eddy simulation

2018

... 基于Q准则的NACA23012-QR冰脊翼型瞬时流场特征(Re_c=2.0×10⁶,Ma=0.21,α=5° ),如图4所示. 其中:Q=0.01,为流场中速度梯度张量的第2矩阵不变;涡结构采用马赫数进行着色.从图4中可以清晰地看出冰脊分离流的空间演化特征,相比于传统的角状冰^[12],冰脊在前缘提前诱导出流动分离,在冰脊前方产生驻涡并诱导出大速度的剪切层流动.这种流动呈现出剪切层Kelvin-Helmholtz不稳定性演化出的管状发卡涡结构,向空间各向扭曲,最终演化为湍流漩涡,与角状冰的涡脱落和破碎现象有一定相似之处.由于前方驻涡的影响,冰脊流的滞止线与翼型切面的夹角增大,使得分离流相比角冰更难以完全再附到翼型表面,这很可能是其在较小攻角下发生尾流剧烈不稳定振荡的主要原因. ...

... t₀,t₁,t₂ 时刻的瞬时压力脉动云图如图10所示.每个时刻间隔1个时间步长,从t₀~t₂分离涡逐渐向下游移动,前序分离涡在发展中逐渐变大.图11提取了P1,P2,P3,P4这4个位置的压力脉动能谱.P1点的压力脉动能谱存在St=0.1325及St=0.025 两个峰值;文献[12]在研究角冰时将其归因于剪切层的涡配对和脱落,这在本研究中同样成立.P2和P3点位置只存在St=0.025的峰值,说明此时剪切层已演化为涡脱落;P4点位于剪切层尾迹中,只能检测到宽频信号,此时剪切层已经演化为低能态大尺度结构. ...

基于IDDES方法的翼型结冰失速分离流动数值模拟

2016

基于IDDES方法的翼型结冰失速分离流动数值模拟

2016

Study on vortex shedding mode on the wake of horn/ridge ice contamination under high-Reynolds conditions

2019

Unsteady modes in flowfield about airfoil with horn-ice shape

2015

... 在最新的试验研究中,文献[15]发现冰脊流场中存在一种区别于常规模态的低频模态,表现为剪切层挥舞形态,与压力波动存在关联.剪切区的谱分析结果表明,这种振荡特性的Strouhal数在0.02附近^[16],与低雷诺数下后台阶流的流动特性比较相似.结合后台阶流理论^[17]不难发现,剪切层尾迹中的涡破裂,引发流体向外卷吸速度和向内的回注速度同时发生变化,这种“流量夹带”和“卷吸流量回注”之间的瞬时不平衡造成了剪切层发展的不稳定,剪切层在这种作用下会伴随尾流发生低频自维持振荡,然而这种现象的本质并没有完全得到解释. ...

Characterization of low-frequency oscillations in the flowfield about an iced airfoil

2014

Analysis of flapping motion of reattaching shear layer behind a two-dimensional backward-facing step

2017

IDDES simulation of flow separation on an 3-D NACA23012 airfoil with spanwise ridge ice

2018

... 本文的研究目的在于采用IDDES方法^[18],对溢流冰脊诱导剪切层振荡的非定常流特性进行高分辨率模拟,期望通过结合正交分解方法,进一步获取近失速形态下,造成冰脊分离流压力波动的主导模态,深入探讨冰脊分离流动不稳定性的时域和频域特征,探索挖掘大尺度相干结构与升力波动的关联性.本文的研究揭示了近失速形态下冰脊分离流导致升力剧烈振荡的物理机制. ...

High-resolution multi-code implementation of unsteady Navier-Stokes flow solver based on paralleled overset adaptive mesh refinement and high-order low-dissipation hybrid schemes

2017

... 计算模型的几何外形如图1所示,其中:

y / c 为 y

方向位置与弦长之比.采用基于CFL3D框架开发的高分辨率流场数值求解器^[19]实现三维冰脊流场的计算模拟.基于连续性方程、动量方程和理想气体方程实现ρ、p、u、v、w、T这6个未知量的计算封闭,其中:ρ为密度;p为压力;u,v,w分别为x,y,z方向的速度;T为温度.忽略体积力和外部热源,其守恒形式可写成矢量形式如下: ...

... 数值离散方面,空间上对流项计算采用基于5阶WENO和6阶中心差分的混合格式进行离散^[19],黏性项采用4阶中心格式;时间上采用双时间步进策略^[20]和隐式近似因子分解内迭代法以2阶精度推进求解.物理时间导数项采用2阶精度的三点向后差分格式离散,伪时间导数项采用1阶向后差分格式逼近. ...

Time dependent calculations using multigrid, with applications to unsteady flows past airfoils and wings

1991

Modal analysis of 3-D iced-airfoil aerodynamics based on proper orthogonal decomposition techniques

2019

... 为了进一步分析冰脊剪切流的不稳定性特征,采用快照正交分解(POD)方法^[21]对流场的压力脉动时域信息进行处理.POD的基本思想是将时间和空间上连续的物理量分解.给定M个时刻的流场快照,任意脉动量可以写作M维POD基和对应模态系数的乘积,则有: ...

〈

〉