力学学报  2018 , 50 (6): 1356-1367 https://doi.org/10.6052/0459-1879-18-129

高速空气动力学

各向同性湍流通过正激波的演化特征研究1)

洪正2), 叶正寅

西北工业大学航空学院,西安 710072

STUDY ON EVOLUTION CHARACTERISTICS OF ISOTROPIC TURBULENCE PASSING THROUGH A NORMAL SHOCK WAVE1)

Hong Zheng2), Ye Zhengyin

School of Aeronautics, Northwestern Polytechnical University, Xi'an 710072, China

中图分类号:  O322,V215.3

文献标识码:  A

收稿日期: 2018-04-18

接受日期:  2018-04-18

网络出版日期:  2018-11-18

版权声明:  2018 力学学报期刊社 力学学报期刊社 所有

基金资助:  1) 国家自然科学基金资助项目(11732013).

作者简介:

2) 洪正,硕士研究生,主要研究方向:激波湍流相互干扰. E-mail: hongzheng@mail.nwpu.edu.cn

展开

摘要

激波与湍流相互作用(shock-turbulence interaction,STI)是空气动力学研究中的一个基础问题.基于格心有限差分法(cell-centered finite difference method,CCFDM)求解器Helios,采用五阶加权紧致非线性格式(weighted compact nonlinear scheme,WCNS)对各向同性湍流通过正激波的情形进行直接数值模拟(direct numerical simulation,DNS).对湍流相关物理量进行统计,分析结果表明,在湍流中波后的密度、温度和压力较无湍流情形下略小,而速度则略大,均在波后呈现短暂过冲然后缓慢向理论值逼近的变化趋势;波后流向雷诺应力突降随之快速增长又衰减,呈现非单调变化趋势,线性相互作用分析(linear interaction analysis,LIA)将其归结为波后能量从声模式转移为涡模式方式,与流向不同,横向雷诺应力突增后单调衰减,波后雷诺应力各向异性明显且随下游距离逐渐增强;波后湍动能突增后呈现非单调变化趋势;泰勒微尺度和Kolmogorov尺度过激波后均明显减小,说明波后湍流长度尺度变小,从而对波后网格的分辨率提出了更高的要求;密度、温度和压力过激波后脉动均方根均增加,密度和压力脉动强度减小,温度脉动强度增大.

关键词: 激波湍流相互作用 ; 格心有限差分法 ; 直接数值模拟 ; 统计分析

Abstract

Shock-turbulence interaction is a kind of important fundamental problem in aerodynamics. Based on solver Helios which applies cell-centered finite difference method (CCFDM), using fifth-order weighted compact nonlinear scheme (WCNS), we conducted direct numerical simulation (DNS) of the situation where isotropic turbulence passes through a normal shock wave. Turbulence statistics are calculated for analysis. We found after shock, density is a little lower than its non-turbulent value, so do temperature and pressure, on the contrary, longitudinal velocity is a little higher than its non-turbulent value. The commonality is that they all show an overshoot immediately behind the shock, after that they gradually approach towards their non-turbulent values along with downstream distance. Longitudinal Reynolds stress suffers a sudden decrease and increases rapidly followed by decaying. This evolution characteristics is captured in linear interaction analysis (LIA) and a transfer of energy from acoustical to vertical modes behind the shock is thought to be accounted for it according to this analysis. Different from longitudinal Reynolds stress, Transverse Reynolds stress suffers a sudden increase then decay monotonically. Anisotropy of Reynolds stress is apparent after shock, and it gradually increases as downstream distance increases. Turbulent kinetic energy suddenly increases and then evolves non-monotonically. Taylor microscale and Kolmogorov scales apparently decrease after shock, indicating the decrease of turbulent length scales, which leads to a requirement of higher resolution of mesh in this zone to solve the flow field. After shock, the root-mean-squares of density, temperature and pressure fluctuations are enhanced, and intensities of density and pressure decrease while intensity of temperature increases.

Keywords: shock-turbulence interaction ; cell-centered finite difference method ; direct numerical simulation ; statistics analysis

0

PDF (6593KB) 元数据 多维度评价 相关文章 收藏文章

本文引用格式 导出 EndNote Ris Bibtex

洪正, 叶正寅. 各向同性湍流通过正激波的演化特征研究1)[J]. 力学学报, 2018, 50(6): 1356-1367 https://doi.org/10.6052/0459-1879-18-129

Hong Zheng, Ye Zhengyin. STUDY ON EVOLUTION CHARACTERISTICS OF ISOTROPIC TURBULENCE PASSING THROUGH A NORMAL SHOCK WAVE1)[J]. Acta Mechanica Sinica, 2018, 50(6): 1356-1367 https://doi.org/10.6052/0459-1879-18-129

引 言

激波与湍流相互作用(shock-turbulence interaction,STI)是一类重要的基础问题,深入理解其中蕴含的物理机理,对于超声速和高超声速转捩、冲压发动机内部流场的演化、燃料的混合以及超声速燃烧过程的理解具有重要意义.在这些应用中,激波与湍流的作用是相互的,它们之间的耦合非常强,包含了复杂的线性和非线性作用机制,引起湍流的结构和统计特性以及激波特性的显著变化.由于涉及的长度、时间尺度范围很广,以及激波和湍流问题中存在的复杂的其他效应,使得这项工作的研究困难重重.

Kovansznay发现对于可压缩湍流中的任意脉动可以分解为3种独立的线性波模型,分别为声波、熵波、涡波[1].基于Kovansznay的分解,Rinber[2-4]使用线性分析的方法对STI问题进行了理论分析,发展出描述该问题的理论模型,称之为线性相互作用分析(linear interaction analysis,LIA),LIA假设激波前湍流脉动相比于平均流动是个小扰动,所以非线性和黏性效应就可以从分析中去除,这样就可以单独分析不同波通过激波后的变化.LIA成功地预测了通过激波后湍动能的增大和湍流长度尺度的减小.Lele[5]利用快速变形理论(rapid distortion theory,RDT)和气体动力学理论推导了在湍流中的激波跳跃关系式,分析表明在特定的静压升下,湍流中的激波速度比在层流中稍大,压缩性增强.

早期Lee等[6-7]和Mahesh等[8]开展的关于STI的DNS研究最具有代表性.Lee等[6-7]考虑了激波强度对激波和各向同性湍流相互作用的影响,发现波后的湍动能和横向涡分量均得到增强,波后湍流尺度变小,如果上游湍流强度较强,则会引起激波面的变形,当湍流马赫数$M_{\rm t}> 0.1(M^2 - 1)$时,激波强度的空间分布就会发生显著变化.Mahesh等[8]研究了上游涡量脉动和熵脉动对激波与各向同性湍流相互作用问题的影响,结果表明,如果上游速度脉动和温度脉动呈负相关,则波后的湍动能、涡量和热力学量脉动的增大会得到增强;如果是正相关,则会被抑制.Mahesh基于体积压缩项和斜压项的相对效应给出了这一现象的解释:体积压缩项和斜压项决定了波后涡量的演化,而上游速度脉动和温度脉动的相关性决定了这两项之间是增强还是抑制的关系,进而影响波后相关量的演化.

Larsson和Lele[9-10]的DNS结果发现过激波后湍流长度尺度,如泰勒微尺度和Kolmogorov尺度会变小,这就要求激波后的网格分辨率较波前要更高以保证正确求解黏性耗散.前人的工作存在网格分辨率不够的问题,这导致结果会出现一些定性的差异,例如波后流向涡量的快速增强以及雷诺应力明显的各向异性.由于激波强度和湍流马赫数的不同,激波在湍流中呈现出两种形态:褶皱激波状态,湍流脉动的分布存在间断;破碎激波状态,湍流脉动的分布是光滑的.Larsson和Bermejo-Moreno[11]在保证网格分辨率足够的前提下,研究了流动雷诺数,平均流动和湍流马赫数对激波与各向同性湍流相互作用的影响.结果表明对于湍动能过激波后的增加,DNS结果与线性理论符合良好;波后湍流各向异性的计算结果与线性理论有定性上的差异;湍流中密度和压力过激波的增大比不考虑湍流的Rankine-Hugoniot值要小并且减小的量与湍流强度的平方有关;特别研究了强湍流下的破碎激波状态,该状态下,局部激波可用光滑压缩替代,激波跳跃概率密度的网格收敛性结果证明这种效应是物理的,并不是数值产生的假象.

Ryu和Livescu[12]在湍流马赫数$1.1 < M_{\rm s} <2.2$和泰勒雷诺数$10 < Re_\lambda <45$参数范围内进行了一系列数值试验.结果表明,Kolmogorov尺度与激波厚度的比值越大,计算结果越接近LIA结果,而比值增大是通过减小湍流马赫数实现的,说明LIA在低湍流马赫数下是可信的,即使上游雷诺数较低的情况.Livescu和Ryu[13]尝试通过LIA代替直接求解激波以减小计算量,从而进行任意平均流马赫数下波后的研究,但是远离激波或非线性及黏性作用明显的区域,未出现期望的衰减.Quadros等[14-15]利用LIA和DNS数据具体研究了波后湍流能量通量,基于线性理论结果提出了新的限制器模型,发展了基于物理的湍动能输运方程模型,对比发现在一系列马赫数范围下模型预测结果与已有DNS结果符合良好.

Tian等[16-17]利用LIA跟DNS研究比较上游多组分流动跟单组分流动在这类问题中的差异,以此来分析密度变化在该问题中产生的效应.发现上游较强的密度变化会使得湍流过激波后的增强及湍流长度尺度的减小相比单组分流动更加显著,激波引起的湍流混合增强也得到了提升,波后混合非对称性更强.

Gao等[18]将五阶WENO(weighted essential non-oscillatory)非线性格式和六阶中心差分线性格式混合,研究了标量场在不同施密特数下随各向同性湍流通过正激波后统计量的演化,并与时间衰减各向同性湍流中的情形进行对比.Vemula和Sinha[19]发现现有的雷诺应力模型在经典的激波与各向同性湍流相互作用问题中预测误差很大,于是他们从物理的角度出发修正了雷诺应力输运方程,结果表明在一系列马赫数下与DNS结果较为符合.Boukharfane等[20]使用高分辨率格式对不同来流马赫数下通过激波后湍流中被动标量的混合进行数值模拟,发现激波显著改变了被动标量场的拓扑结构,增强了标量混合,且随来流马赫数增大而愈发明显.

王国蕾和陆夕云[21]对激波与湍流相互作用问题研究进展进行了总结和展望,目前国内关于STI问题的研究主要集中在湍流边界层与激波相互作用方面[22-24],对于各向同性湍流与激波相互作用问题的研究几乎是空白.

本文基于课题组近期发展的高阶精度数值计算方法和程序[25-26],使用DNS方法对各向同性湍流与激波相互作用问题进行数值模拟,统计流场信息进行分析,观察湍流统计量过激波后的变化并进行分析.

1 控制方程及数值方法

1.1 控制方程

DNS方法的实质是求解非定常的Navier-Stokes方程,不使用任何模型,直接对流动中的各种大小尺度流动进行求解,要求网格尺度与时间步长要与流动最小尺度相当,相比于雷诺平均(Reynolds-averaged,RANS)方法和大涡模拟(large eddy simulation,LES)方法,对计算资源提出了更高的要求.本文使用无量纲的Navier-Stokes方程以及理想气体方程

$$\frac{\partial \rho }{\partial t} + \frac{\partial }{\partial x_i}(\rho u_i ) = 0 (1)$$$$\frac{\partial }{\partial t}(\rho u_i ) + \frac{\partial}{\partial x_j }(\rho u_i u_j + p\delta _{ij} ) - \frac{\partial\tau _{ij} }{\partial x_j } = 0(2)$$$$\frac{\partial }{\partial t}(\rho E) + \frac{\partial }{\partial x_i }((\rho E + p)u_i ) - \frac{\partial }{\partial x_j }(u_i \tau_{ij} ) + \frac{\partial q_i }{\partial x_i } = 0(3)$$$$\frac{1}{\gamma }\rho T (4)$$

其中

$$\tau _{ij} = \frac{\mu }{Re}\left(\frac{\partial u_i }{\partial x_j } + \frac{\partial u_j }{\partial x_i } -\frac{2}{3}\frac{\partial u_k }{\partial x_k }\right) \\ q_i = - \frac{\mu c_p }{(\gamma - 1)RePr }\frac{\partial T}{\partial x_i } \\ E = \frac{p}{\rho (\gamma - 1)} + \frac{1}{2}u_i u_i \\ \gamma = \frac{c_p }{c_v } = 1.4$$

式中,$ \mu = \left(\dfrac{1 + C}{T +C}\right)T^{1.5},C = 110.4K / T_0 ,T_0 = 273.15{\rm K}$

无量纲特征参数

$$Re = \frac{\rho _0 c_0 L}{\mu _0 }$$

$$Pr = \frac{\mu _0 c_p }{k_0 }$$

1.2 数值格式

本文基于课题组开发的格心有限差分求解器Helios,在满足几何守恒率的条件下实现高阶精度格式,相比于传统的格点有限差分求解器,格心有限差分求解器对于多块网格边界的处理更简单,降到二阶后等价于有限体积法,通量分裂格式兼容矢通量分裂格式和通量差分裂格式[25-26].

本文使用邓小刚提出的五阶加权紧致非线性格式(weighted compact nonlinear scheme,WCNS)计算无黏项,该格式在光滑区域近似为五阶紧致迎风格式,在激波附近降阶抑制非物理的振荡[27],黏性项采用六阶紧致中心差分格式计算,通量分裂使用Roe格式. 时间推进采用经典的四阶Runge-Kutta显式格式.

1.3 入口条件

平均流动参数为$\overline {u_1 } = U_1,~\overline {u_2 } =\overline {u_3 } = 0$, $\overline p =1/\gamma $, $\overline \rho =1$.对于脉动信息,首先需要计算一个时间衰减各向同性湍流场,速度导数偏斜因子随着时间从零开始逐渐减小,通常认为速度导数偏斜因子在$-0.4\sim-0.6$之间表明各向同性湍流达到充分发展状态[29],根据计算需求选择充分发展状态的某一时刻流场作为入口湍流数据源.

具体做法是:从选定时刻的流场提取出速度、压力和密度脉动值,然后将物理空间中的脉动信息变换到傅里叶空间,离散傅里叶正变换表示如下

\begin{equation}\label{eq5} \mathop f\limits^ \wedge (k_1 ,k_2 ,k_3 ) =\frac{1}{N_x N_y N_z }\sum\limits_{x,y,z} {f(x,y,z){\rm e}^{ -{\rm i}(k_1 x + k_2 y + k_3 z)}} \tag{5}\end{equation}

其中,$k_1, k_2, k_3$是谱空间波数矢量,$N_x ,N_y ,N_z$是三个方向离散网格数,$x,y,z$是物理空间中的坐标.

利用泰勒假设冻结流场,湍流脉动以平均对流速度通过边界,即通过$x =x_0 + U_1 t$,不失一般性,$x_0$取为零,将来流方向$x$轴转换为入口上的时间$t$轴,离散傅里叶反变换为

$$ f(t,y,z) = \sum\limits_{k_1 ,k_2 ,k_3 } {\mathop f\limits^ \wedge } (k_1 ,k_2 ,k_3 ){\rm e}^{{\rm i}(k_1 U_1 t +k_2 y + k_3 z)} (1)$$

前人的工作表明,湍流马赫数较低的情况下该假设是有效的.因为是超声速入口,不同时刻对应的速度、压力和密度脉动值与平均流动叠加之后直接在入口给定即可.

1.4 出口条件

激波后是亚声速流动,避免出口边界产生非物理的压力反射是需要慎重考虑的.通过增加吸收层,将湍流出口逐渐过渡到层流出口,以此来消除影响[9].具体做法是在吸收层的控制方程右侧增加源项:$ - \sigma ({(x -x_{sp}) }/{(x_{\max } - x_{sp} )})^2(f - \overline {f_{yz} })$其中,$x_{sp} $和$x_{\max }$分别是吸收层的起始和终止$x$轴位置,$f = \rho ,~\rho u_i ,~\rho e_0 $表示守恒变量,$\overline {f_{yz} }$表示在$yz$平面上的平均,$\sigma = {k_0 u_1 }/{2\pi}$,是个常数,本文中$x_{sp} = 3\pi $, $x_{\max } - x_{sp} = \pi$.

2 数值方法验证

首先需要验证施加湍流入口条件所需单独计算的入口湍流场即时间衰减各向同性湍流所使用数值方法的合理性.其次得到入口湍流场后,利用泰勒假设设置湍流入口的有效性需要验证,与时间衰减湍流对应,称之为空间衰减各向同性湍流.最后需要验证五阶WCNS格式对于激波的捕捉能力以及精度.

2.1 时间衰减各向同性湍流

考虑计算域为笛卡尔坐标系下的$[0,2\pi ]\times [0,2\pi ]\times[0,2\pi]$立方体,该区域内充满可压缩的黏性流体,初始时刻流体存在复杂的运动,在没有外力的作用下,该区域内的流动会在黏性的作用下逐渐衰减.3个方向均匀各向同性,均采用周期性边界条件.初始流场速度需要满足给定的一维能谱$E(k) = Ak^4\exp ( - 2k^2 /k_0^2 )$,其中$k$是波数,$k_0$是能谱峰值波数,$A$是与初始动能相关的常数,同时满足不可压质量守恒方程即散度为零的条件,使用Rogallo提出的方法初始化速度场[30-31],压力和密度设置为常值,平均流动速度为零.

相关参数定义如下:泰勒微尺度$\lambda ^2 ={u'^2}/$ ${\overline {(\partial u_1 / \partial x_1 )^2}}$,基于泰勒微尺度的雷诺数$Re_\lambda = {u'\lambda \overline \rho}/{\overline \mu }$,湍流马赫数$M_t = {\sqrt 3 u'}/{\overline c}$,湍流脉动速度$u'^2 = \overline {(u_1^2 + u_2^2 + u_3^2 )} $.计算参数选择如下:$A = 0.00013$, $k_0 = 8$, $Re_\lambda = 72$,$M_t = 0.5$, $Pr = 0.72$. 时间步长$\Delta t = 0.001$.网格为$128^3$均匀网格.

图1为时间发展过程中归一化湍动能$K(t) / K_0$及速度导数偏斜因子的变化情况,其中湍动能的定义为:$K(t) =\overline {\rho u_i u_i } / 2$,$K_0$为初始时刻的湍动能,速度导数偏斜因子定义为$S_3 ({\partial u_1}/{\partial x_1 }) = {\overline {(\partial u_1/\partial x_1 )^3} }/{\overline {(\partial u_1 / \partial x_1 )^2} ^{3 / 2}}$.图中横坐标为无量纲时间$t / \tau $,$\tau = \sqrt{({32}/{A})^{1/2}}/(2\pi )^{1 / 4}K_0^{ - 7 / 2}$为初始时刻的大涡转换时间(large-eddy-turnover time).从图1可以看出,本节计算的湍动能衰减情况与Samtaney等[32]的结果吻合十分理想,此外,速度导数偏斜因子与文献符合良好,偏斜因子是一个高阶统计量,对流场中的小尺度流动十分敏感,说明该网格分辨率下,迎风格式带来的数值耗散和物理耗散相比是可以忽略的.

图1   平均湍动能及速度导数偏斜因子随时间的变化

Fig.1   Evolution of averaged turbulent kinetic energy and velocity derivative skewness

图2为湍流马赫数和归一化的扰动密度均方根随无量纲时间的变化.\vspace{2mm}归一化的密度脉动均方根定义为$\rho _{\rm rms} = \sqrt{\overline {\rho '^2} } / M_t^2 (0)$,其中$\rho ' = \rho -\overline \rho $为密度脉动.与文献中的结果对比可以看出,这两个量与文献吻合得十分理想.

图2   湍流马赫数及归一化的脉动密度均方根随时间的变化

Fig.2   Evolution of the turbulent Mach number and normalized root-mean-square density fluctuation

图1图2说明了本文采用的有限差分求解器Helios和五阶WCNS格式对于本湍流的直接数值模拟是有效的.

2.2 空间衰减各向同性湍流

与时间衰减各向同性湍流相比,空间衰减各向同性湍流的难点是在于流向的不均匀性,在这个方向上,因为入流为超声速,利用泰勒假设将脉动直接施加在流动入口边界处,每一时间步更新入口值.湍流数据来自单独计算的某一时刻的时间衰减各向同性湍流场,该时刻湍流场为充分发展状态.出流为超声速,出口值由内场插值得到.计算域为笛卡尔坐标系下$[0,2\pi ]\times [0,2\pi ]\times[0,2\pi ]$立方体,网格为$128^3$的均匀网格. 计算参数:$Re \approx68$, $Pr = 0.72$,$M_1 = 1.5$, $M_t \approx 0.104$, $Re_\lambda\approx 5.6$,时间步长$\Delta t = 0.001$.

时间各向同性湍流随时间衰减,空间各向同性湍流随空间距离衰减,为了能对两者进行比较,利用$t= x_1 / U_1$将空间位置转换为时间,用时间湍流初场的大涡转换时间$\tau$进行无量纲化.图3所示为时间/空间衰减各向同性湍流湍流马赫数随无量纲时间的变化,可以看到两者符合良好.速度导数偏斜因子$S_3$是判断各向同性湍流是否真实且处于完全发展状态的重要依据,图4表明时间/空间衰减各向同性湍流速度导数偏斜因子均在$-0.4\sim-0.6$区间内且较为一致.上述结果表明在本节的计算条件下,利用泰勒假设设置湍流入口边界是有效的至于五阶WCNS格式对激波的捕捉及精度的验证,在文献[24,25]中有详细的算例验证过该格式的有效性,此处不再赘述.

图3   湍流马赫数随衰减曲线

Fig.3   Evolution of the turbulent Mach number

图4   速度导数偏斜因子变化曲线

Fig.4   Evolution of the velocity derivative skewness

2.3 网格收敛性验证

湍流的DNS计算需要完全求解粘性耗散,并且在激波后湍流的Kolmogorov尺度会减小,故波后网格的分辨率比波前要求更高,结果对于网格的依赖性需要进行验证.相同参数下,在一系列逐步加密的网格上进行计算,流向和横向的雷诺应力随空间位置变化的结果分别如图5图6所示.

图5   流向雷诺应力随空间位置的变化

Fig.5   Evolution of the longitudinal Reynolds stress

图6   横向雷诺应力随空间位置的变化

Fig.6   Evolution of the transverse Reynolds stress

图例中$108\times 32\times32$表示流向采用108个网格点,剩余两个方向采用32个网格点,其他网格以此类推.从结果来看,除了$108\times 32\times32$网格结果与其他网格结果有明显偏差外,其余网格结果基本一致,说明本文结果使用的$302\times128\times 128$网格的分辨率已经足够描述流场.

3 DNS结果及分析

3.1 计算模型

图7为计算模型示意图,计算域为$(L_1 ,L_2 ,L_3 ) = (4\pi ,2\pi,2\pi)$,$x$轴正方向为平均来流方向,$y,z$轴方向均匀各向同性,采用周期性边界条件.$x = 0$位置为超声速入口,湍流信息在入口给定;$x = \pi$位置处利用R-H关系初始化得到一道正激波;$x = 3\pi $到$x = 4\pi$处是吸收层,将湍流出口过渡到层流出口,以消除出口处的非物理压力反射带来的误差.流向即$x$轴方向上在激波处局部加密网格,以保证激波厚度内至少有10个网格,吸收层采用拉伸网格,一方面增强湍流的衰减,另一方面减少计算量.网格如图8所示,$y,z$轴方向使用均匀网格,网格量为$302\times128\times 128$. 入口参数为$Re \approx 68$, $Pr = 0.72$,$M_1 =1.5$, $M_t \approx 0.104$, $Re_\lambda \approx 5.6$.时间步长$\Delta t = 0.001$.

图7   计算模型示意图

Fig.7   The sketch of computational model

图8   计算网格

Fig.8   Computational mesh

当入口处最初的湍流脉动发展到出口处时,过渡阶段结束,在接下来的$2\tau$时间内取60个样本流场,在各向同性的$y,z$轴方向进行平均计算统计量.流向即$x$轴距离有多种方式可以进行无量纲化,本节选择能谱峰值波数$k_0$对流向距离进行无量纲化. 所有结果图片中,平均激波位置被平移到$x =0$处,这种处理之下,入口处的位置为$k_0 x \approx -25$,出口位置为$k_0 x \approx 75$,吸收层起始位置为$k_0 x \approx50$,$y,z$轴方向长度$k_0 x = 16\pi $.统计量的平均在均匀的$y,z$轴方向进行的,雷诺平均表示为$\overline f$,密度平均表示为$\widetilde{f} = {\overline {\rho f} }/{\overline\rho }$,相应的脉动分量分别表示为$f' = f - \overline f $和$f" = f- \widetilde{f}$.本节图中灰色区域为激波非定常运动区域,该区域平均散度为负值,表现为压缩效应.

3.2 平均流场

图9$\sim$图11分别为平均密度、温度和压力剖面,三者变化趋势类似,过激波后先跳升,紧接着少量减小,然后随着下游距离缓慢向由R-H关系得到的值逼近,过程中一直低于理想值.图12为流场平均流向速度剖面,过激波后先陡降,紧接着少量增大,然后随着下游距离缓慢向由R-H关系得到的值逼近,过程中一直高于理想值.

图9   平均流场密度剖面图(长虚线为R-H值)

Fig.9   Profile of the mean density (long-dash line denotes R-H value)

图10   平均流场温度剖面图

Fig.10   Profile of the mean temperature

图11   平均流场压力剖面图

Fig.11   Profile of the mean pressure

图12   平均流场流向速度剖面图

Fig.12   Profile of the mean longitudinal velocity

平均来说,流动被激波压缩后,紧跟着有轻微的膨胀.

Lele[1]推导了湍流中的激波跳跃关系式,得到了湍流中激波过激波后压力增加较无湍流情况下变小的结论,这与数值的结果是一致的.

3.3 雷诺应力

图13图14分别显示了过激波前后流向和横向雷诺应力的变化情况,用激波前的相应值进行无量纲化.图中可以看出,流向雷诺应力过激波后突然减小到0.7左右,然后快速增加到峰值1.2附近,随后单调衰减;横向雷诺应力通过激波后突然增强到1.5附近,随后单调衰减,LIA也捕捉到了这一非单调的趋势,并把原因归结为波后能量从声模式转移到涡模式[10].图15中湍动能在通过激波后得到了增强,并呈现出非单调的变化,这也是流向和横向雷诺应力波后变化的综合效果.需要说明的是,Larsson 等[11]通过简单的数值模型指出平均激波位置内的峰值主要是由于激波的非定常运动引起的.

图13   流向雷诺应力随空间位置的变化

Fig.13   Evolution of the longitudinal Reynolds stress

图14   横向雷诺应力随空间位置的变化

Fig.14   Evolution of the transverse Reynolds stress

图15   湍动能随空间位置的变化

Fig.15   Evolution of the turbulent kinetic energy

图16从流向雷诺应力与横向雷诺应力的比值来看,激波前比值在1附近,流向和横向是各向同性的,但是穿过激波之后,比值先降到0.5附近,然后快速上升到1.2左右,而后随着下游距离的增加以较小的速度继续增加.这反映经过激波后,雷诺应力的各向同性被破坏了,流向和横向脉动出现了明显的差异,且随着下游距离各向异性逐渐增强.

图16   流向与横向雷诺应力之比随空间位置的变化

Fig.16   Evolution of the ratio of longitudinal Reynolds stress to transverse one

3.4 涡量

涡量通过激波的变化情况见图17$\sim$图19,用激波前的相应值进行无量纲化.流向涡量通过激波后没有明显的增强或者减弱,只是在激波内部出现峰值随即下降;而横向涡量通过激波后突增紧接着单调衰减.平均来说,流向涡量平行于激波面,几乎不受压缩,故通过激波没有明显变化,而横向涡量垂直于激波面,受到明显的压缩,从流向和横向涡量分量比值来看,激波前各向同性,比值约为1,过激波后,横向涡量增强为流向的3倍左右,随着下游位置先快速减小后缓慢减小,到达吸收层前约为1.5.这说明涡量在通过激波后,涡量是各向异性的,从趋势看各向异性随下游距离逐渐减弱,但是没有消除.

图17   流向涡量随空间位置的变化

Fig.17   Evolution of the longitudinal vorticity

图18   横向涡量随空间位置的变化

Fig.18   Evolution of the transverse vorticity

图19   横向与流向涡量比值随空间位置的变化

Fig.19   Evolution of the ratio of longitudinal vorticity to transverse one

3.5 湍流长度尺度

图20反映了通过激波前后泰勒微尺度的变化情况,流向和横向尺度的定义分别为

激波前,流向和横向泰勒微尺度大致保持一致,随着湍流的发展逐渐增加;激波后,流向和横向尺度均有明显减小,但流向减小得更加剧烈,接着流向尺度先快速增长然后保持一个较缓的速率持续增长;横向尺度先快速增大达到峰值然后稍有减小,而后随下游距离均匀增大.穿过激波后,横向泰勒微尺度的增长速率没有明显变化,而流向尺度的增长速率先增大后减小,减小后的增长速率仍比波前大.

图20   流向和横向泰勒微尺度随空间位置的变化

Fig.20   Evolution of the longitudinal and transverse Taylor microscale

LIA分析得到的结果表明通过激波后$\lambda _1 < \lambda _2$,Lee等[7]的DNS结果也呈现了相同的结论,但是Larsson网格分辨率、格式精度更高的DNS结果却表明,激波后$\lambda_1 < \lambda _2 $,但是$\lambda _1 $迅速增加以致短距离后$\lambda_1 > \lambda _2 $.Larsson认为,与LIA结果相反的原因是由于在激波后泰勒微尺度的演化过程中非线性作用至关重要,而LIA忽略了非线性影响,Lee的结果由于网格分辨率不够,非线性的预测不足,所以更倾向于线性的结果[10].本文的结果是激波后距离约20内$\lambda _1 < \lambda _2$,之后$\lambda _1 > \lambda _2 $,趋向于Larsson的结果.

图21中Kolmogorov尺度定义为$\eta = ((\mu / \rho )^3 / \varepsilon)^{1 / 4}$, $\varepsilon $为湍动能耗散率.该尺度被认为是流场中最小的尺度,湍动能在该尺度下通过黏性作用被转化为流体内能.Kolmogorov尺度通过激波后减小了约1/3,无间断时均随距离线性增长.DNS要求直接求解粘性耗散,故对计算网格量分辨率有一定的要求,文献${[33]}$指出一般需要满足$k_{\max} \eta \ge 1.5$, $k_{\max }$为可分辨的最大波数,为单方向网格量的一半.对于本文的计算状态,$\eta _{\min } \approx 0.05$, $k_{\max } =64$, $k_{\max } \eta _{\min } \approx3$足够满足网格分辨率的要求,故可以说本文的计算结果是合理且可信的.

图21   Kolmogorov尺度随空间位置的变化

Fig.21   Evolution of the Kolmogorov scale

3.6 热力学量

图22为密度、温度和压力脉动均方根随空间位置的变化,其中脉动均方根定义为$f_{\rm rms} = (\overline {(f - \overline f )^2} )^{1 / 2}$, $f = \rho,T,P$. 用对应的波前值进行无量纲.对比可以看出,密度和温度脉动的变化趋势十分相似,波前随空间距离湍流逐步衰减,密度和温度脉动随之逐渐减小,通过激波后,脉动均方根突增后快速衰减到2附近之后随下游距离缓慢衰减;从波后增强效果来看,压力最大,温度最小.从图23中的脉动均方根与平均值的比值也即脉动强度来看,脉动强度过激波后的趋势与均方根一致,但是因为过激波后平均密度、温度和压力会升高,所以幅度变小.总体来看,波后密度和压力的脉动强度减小了,其中压力最为明显,只有温度脉动强度增加了;波前各量的脉动强度相差很大,波后却趋于一致

图22   脉动均方根随空间位置的变化

Fig.22   Evolution of the root-mean-square fluctuation

图23   脉动强度随空间位置的变化

Fig. 23   Evolution of the fluctuation intensity

3.7 流场

图24使用Q准则识别流场中的涡结构,深色面表示瞬时激波面,用脉动压力着色.图中波前涡结构杂乱无章,各方向均匀,而波后涡结构明显不再各向同性,呈现出一定的方向性.图25为$y = \pi$位置截面的脉动马赫数等值线图,图中因为等值线密集而显示为黑色线的是瞬时激波面,由于湍流脉动的存在,激波面是一条曲线,并且局部激波的厚度也随脉动的变化而增厚或者变薄.

图24   Q准则识别的流场涡结构(中间深色平面为瞬时激波面,用脉动压力着色)

Fig.24   Vortex structure characterized by Q criterion (dark plane represents the instantaneous shock surface, colored by pressure fluctuation)

图25   $y=\pi$位置截面脉动马赫数等值线图

Fig.25   Contour of fluctuation Mach number at $y=\pi$ plane

4 结 论

采用DNS方法对各向同性湍流与通过正激波的情形进行研究,通过统计与分析得到以下结论:

(1)同样的平均来流流动参数下,激波后的密度、温度和压力比无湍流情形下要低;波后的流向速度比无湍流情形下要高.

(2)过激波后,流向雷诺应力突降后快速达到峰值后衰减,横向雷诺应力突增后单调衰减.从流向与横向应力比值来看,激波后各向异性明显且随下游距离逐渐增强.湍动能波后得到增强,并呈现非单调的变化趋势.

(3)流向涡量在激波内有小幅上升随即衰减,激波前后没有明显变化;横向涡量波后明显的增强,呈现单调衰减的趋势.横向流向涡量的比值表明波后涡量的各向差异达到最大,然后随着下游距离逐渐减弱,但在本文计算中没有恢复到各向同性.

(4)波前各方向泰勒微尺度随空间位置逐渐增加,过激波后,流向横向均减小,流向减小得更多,随即快速增加,然后以较小速率继续随下游增大;横向尺度经过短暂的上升下降后,以与波前相同的速率随下游增加.波后$\lambda _1 < \lambda _2 $,经过一段距离后则$\lambda _1 >\lambda _2 $.Kolmogorov尺度过激波后减少约1/3,故对波后网格的分辨率提出了更高的要求.

(5)密度、温度和压力脉动均方根在穿过激波后均得到了增强;从脉动强度来看,密度和压力脉动强度波后减小,温度脉动强度增大,波前各脉动强度差异明显,波后却趋于一致.

The authors have declared that no competing interests exist.


参考文献

[1] Leslie SG,

Kovasznay. Turbulence in supersonic flow

. Journal of the Aeronautical Sciences, 1953, 20: 657-682

[本文引用: 2]     

[2] Ribner HS.Convection of a pattern of vorticity through a shock wave. NACA TN 2864, 1953

[本文引用: 1]     

[3] Ribner HS.Shock/turbulence interaction and the generation of noise. NACA TN 3288, 1954

[4] Ribner HS.

Spectra of noise and amplified turbulence emanating from shock/turbulence interaction

. AIAA Journal, 1987, 35: 436-442

[本文引用: 1]     

[5] Lele SK.

Shock-jump relations in a turbulent flow

. Physics of Fluids A : Fluid Dynamics, 1992, 4: 2900-2905

[本文引用: 1]     

[6] Lee S, Lele SK, Moin P.

Direct numerical simulation of isotropic turbulence interacting with a weak shock wave

. Journal of Fluid Mechanics, 1993, 251: 533-562

[本文引用: 2]     

[7] Lee S, Lele SK, Moin P.

Interaction of isotropic turbulence with shock waves: effect of shock strength

. Journal of Fluid Mechanics, 1997, 340: 225-247

[本文引用: 3]     

[8] Mahesh K, Lele SK, Moin P.

The influence of entropy fluctuations on the interaction of turbulence with shock wave

. Journal of Fluid Mechanics, 1997, 334: 353-379

[本文引用: 2]     

[9] Lele SK, Larsson J.

Shock-turbulence interaction: what we know and what we can learn from peta-scale simulations

. Journal of Physics: Conference Series, 2009, 180: 012032

[本文引用: 2]     

[10] Larsson J, Lele SK.

Direct numerical simulation of canonical shock/turbulence interaction

. Physics of Fluids, 2009, 21: 126101

[本文引用: 3]     

[11] Larsson J, Bermejo-Moreno I, Lele SK.

Reynolds- and Machnumber effects in canonical shock-turbulence interaction

. Journal of Fluid Mechanics, 2013, 717: 293-321

[本文引用: 2]     

[12] Ryu J, Livescu D.

Turbulence structure behind the shock in canonical shock-vortical turbulence interaction

. Journal of Fluid Mechanics, 2014, 756: R1

[本文引用: 1]     

[13] Livescu D, Ryu J.

Vortitity dynamics after the shock-turbulence interaction

. Shock Waves, 2016, 26: 241-251

[本文引用: 1]     

[14] Quadros R, Sinha K, Larsson J.

Turbulent energy flux generated by shock/homogeneous-turbulence inter-action

. Journal of Fluid Mechanics, 2016, 796: 113-157

[本文引用: 1]     

[15] Quadros R, Sinha K.

Modeling of turbulent energy flux in canonical shock-turbulence interaction

. International Journal of Heat and Fluid Flow, 2016, 61: 626-635

[本文引用: 1]     

[16] Tian YF, Jaberi FA, Livescu D, et al.

Numerical simulation of multi-fluid shock-turbulence interaction

. AIP Conference Proceedings, 2017, 1793:150010

[本文引用: 1]     

[17] Tian YF, Jaberi FA, Li ZR, et al.

Numerical study of variable density turbulence interaction with a normal shock wave

. Journal of Fluid Mechanics, 2017, 829: 551-588

[本文引用: 1]     

[18] Gao XY, Bermejo-Moreno I, Larsson J.

Direct numerical simulation of passive scalar mixing in shock turbulence interaction

//70th Annual Meeting of the APS Division of Fluid Dynamics, Denver, Colorado, 2017

[本文引用: 1]     

[19] Vemula JB, Sinha K.

Reynolds stress models applied to canonical shock-turbulence interaction

. Journal of Turbuulence, 2017, 18: 653-687

[本文引用: 1]     

[20] Boukharfane R, Bouali Z, Mura A.

Evolution of scalar and velocity dynamics in planar shock-turbulence interaction

. Shock Waves, 2018, 5: 1-25

[本文引用: 1]     

[21] 王国蕾, 陆夕云.

激波和湍流相互作用的数值模拟

.力学进展, 2012, 42(3): 274-281

[本文引用: 1]     

(Wang Guolei, Lu Xiyun.

Numerical simulation of shock wave/turbulence interactions

. Advances in Mechanics, 2012, 42(3): 274-281(in Chinese))

[本文引用: 1]     

[22] 崔光耀, 潘翀, 高琪.

沟槽方向对湍流边界层流动结构影响的实验研究

. 力学学报, 2017, 49(6): 1201-1212

[本文引用: 1]     

(Cui Guanyao, Pan Chong, Gao Qi, et al.

Flow structure in the turbulent boundary layer over directional riblets surfaces

. Chinese Journal of Theoretical and Applied Mechanics, 2017, 49(6): 1201-1212(in Chinese))

[本文引用: 1]     

[23] 童福林, 李欣, 于长平.

高超声速激波湍流边界层干扰直接数值模拟研究

. 力学学报, 2018, 50(2): 197-208

(Tong Fulin, Li Xin, Yu Changping, et al.

Direct numerical simulation of hypersonic shock wave and turbulent boundary layer interactions

. Chinese Journal of Theoretical and Applied Mechanics, 2018, 50(2): 197-208(in Chinese))

[24] 高天达, 孙姣, 范赢.

基于PIV技术分析颗粒在湍流边界层中的行为

. 力学学报, DOI: 10.6052/0459-1879-18-211

[本文引用: 2]     

(PIV experimental investigation on the behavior of particle in the turbulent boundary layer. Chinese Journal of Theoretical and Applied Mechanics, DOI: 10.6052/0459-1879-18-211

[本文引用: 2]     

[25] Liao F, Ye ZY, Zhang LX.

Extending geometric conservation law to cell-centered finite difference methods on stationary grids

. Journal of Computational Physics, 2015, 284: 419-433

[本文引用: 3]     

[26] Liao F, Ye ZY.

Extending geometric conservation law to cell-centered finite difference methods on moving and deforming grids

. Journal of Computational Physics, 2015, 303: 212-221

[本文引用: 2]     

[27] Deng XG, Zhang HX.

Developing high-order weigh-ted compact nonlinear schemes

. Journal of Computational Physics, 2000, 165: 22-44

[本文引用: 1]     

[28] 刘昕, 邓小刚, 毛枚良.

高阶精度非线性格式WCNS-E-5在二维流动中的应用研究

.空气动力学报, 2004, 22(2): 206-210

(Liu Xin, Deng Xiaogang, Mao Meiliang, et al.

Weighted compact high-order nonlinear scheme WCNS-E-5 applied to two dimensional flows

. Chinese Journal of Theoretical Applied Mechanics, 2004, 22(2): 206-210(in Chinese))

[29] Mahesh K, Moin P, Lele SK.The interaction of a shock wave with a turbulent shear flow. AFSOR TF- 69, 1996

[本文引用: 1]     

[30] Rogallo RS.Numerical experiments in homogeneous turbulence. NASA, Technical Report 81315, 1981

[本文引用: 1]     

[31] 秦泽聪, 方乐.

一种改进的均匀各向同性湍流初始化方法

. 力学学报, 2016, 48(6): 1319-1325

[本文引用: 1]     

(Qin Zecong, Fang Le.

An improved method for initializing homogeneous isotropic turbulent flows

. Chinese Journal of Theoretical and Applied Mechanics, 2016, 48(6): 1319-1325(in Chinese))

[本文引用: 1]     

[32] Samtaney R, Pullin DI, Kosovic B.

Direct numerical simulation of decaying compressible turbulence and shocklet statistics

. Physics of Fluids, 2001, 13: 1415-1430

[本文引用: 1]     

[33] Pope SB. Turbulent Flows.Cambridge: Cambridge University Press, 2000

[本文引用: 1]     

/