受限撞击流反应器(confined impinging jet reactor, CIJR)因其快速、高效的混合性能和独特的相间传递特性，在聚合物反应注射成型^[1-2]、旁侧进气突扩燃烧室^[3]、固体颗粒干燥^[4]、纳米材料合成^[5]等工业生产中得到了广泛应用。

近年来，国内外众多学者对不同控制参数下CIJR的流动特性和混合特性进行了详细的实验^[6-9]和数值模拟^[10-12]研究。已有研究结果表明：CIJR的湍流流场中存在大量的拟序涡结构，这些结构对流场的流动特性和混合效果有重要影响。Schwertfirm等^[13]在低雷诺数下对CIJR进行直接数值模拟(director numerical simulation, DNS)研究，发现稳定涡是影响主混合管道内流动的主要因素。Wang和Mujumdar^[14]采用标准k-ε湍流模型对受限湍流圆形对置撞击流的流动和混合特性进行数值研究, 使用温度作为示踪剂并定义混合指数(mixing index, MI)反映混合效果，发现对置射流碰撞产生的2个涡旋控制着圆顶区、撞击区及其附近的混合特性。Shi等^[15]采用平面激光诱导荧光(planar laser induced fluorescence, PLIF)技术研究水和甘油溶液在有、无激励的CIJR中的混合特性，使用离析度(intensity of segregation, IOS)反映混合效果。结果表明：在没有激励的情况下，Re=100时在撞击点上方的圆顶区形成的涡旋可以改善涡接触面附近的混合；Re=150时在撞击点下游形成的大尺度涡进一步增强了混合。Zhang等^[16]采用PLIF技术研究了T形撞击流反应器的三维流动特性和混合特性，同样使用IOS表征混合效果。结果表明：在非定常吞噬流中撞击面上周期性地出现涡合并现象，显著改善了混合效果。以上研究通常以整个反应器为对象，使用IOS、MI等参数直观反映混合效果，重点分析涡结构对宏观混合的影响，而揭示能量输运过程和强化流体混合需要进一步深入细致研究各种涡的生成和演化机制。部分学者通过标量分布来分析流场的混合特性^[17-19]，但并未深入研究涡结构对标量输运过程的影响。

在CIJR的大部分区域，湍流流动是非均匀、各向异性的^[20]，需要可靠有效的数值方法对流场进行求解。DNS或大涡模拟(large eddy simulation, LES)以其提供全四维(时间和空间)动力学特性和求解所有重要湍流尺度运动的能力，被视为研究流动机理的可靠手段。然而，Spalart等^[21]认为对于高雷诺数的边界层流动，LES和DNS需要大量的网格，使得计算成本明显高于雷诺平均方法(Reynolds-averaged navier-stokes, RANS)。改进的延迟分离涡模拟(improved delayed detached-eddy simulation, IDDES)方法^[22]能够实现在近壁区采用RANS求解和在远离壁面区域采用LES求解的自动切换，适用于求解高雷诺数下的复杂湍流流动。在过去的几年中，IDDES方法通常被用来研究高速列车尾流动力学^[23-25]和风力机翼型的分离流动^[26]，体现了IDDES方法在尾流动力学建模中的能力。另外，一些研究人员使用分离涡方法(detached-eddy simulation, DES)和延迟分离涡方法(delayed detached-eddy simulation, DDES)模拟T形管的湍流混合和热混合特性^[27]、明渠流动中的拟序结构^[28]和单射流撞击壁面的换热特性^[29-30]。这些研究提高了我们使用IDDES方法模拟CIJR内部流场的信心。

本文使用IDDES方法对CIJR内部流场进行三维数值模拟，采用相平均方法对脉动信号进行分解，通过分析相平均速度场和瞬时速度场揭示流场不同区域涡结构的动力学行为对流动特性的影响；最后利用相平均温度场和热通量场讨论涡结构对湍流混合的影响。

1 模型构建 1.1 物理模型

图 1给出了所研究的CIJR结构示意图，几何模型和相应的尺寸基本与Liou等^[20]的实验装置一致。参考已有的研究工作^{[20, 31-32]}，将计算域划分为4个区域：自由射流区、圆顶区、撞击区和下游发展区。采用空气作为模拟介质，2股射流分别从上下2个同轴、相向对置的矩形喷嘴射入长方体腔室，由于射流喷嘴附近的突扩效应和速度剪切，以及射流撞击引起的扩散，在圆顶区产生一对反向旋转的涡旋，称为一次回流区(primary recirculation zones, PRZ)；2股射流相互撞击后旋转90°向下游发展，在喷嘴的下游边缘角点分离形成2个剪切层，之后分离剪切层在上、下壁面再附着，形成另一对反向旋转的涡旋，称为二次回流区(secondary recirculation zones, SRZ)^{[20, 33]}；最终2股流体从腔室下游侧的出口流出。如图 1，在笛卡尔坐标系中，X、Y、Z分别代表流向、横向和展向坐标，U、V、W为相应的速度分量。流向坐标零点位于喷嘴的上游边缘，其他2个方向的坐标零点分别位于腔室的横向和展向中心。矩形喷嘴的横向截面尺寸为L_S×W_C=15 mm×120 mm，长方体腔室的流向截面尺寸为H_C×W_C=30 mm×120 mm，圆顶高度为H_D=0.5H_C，腔室下游侧的长度为L_C=9H_C。

1.2 数学模型

为简化数学模型，本文引入以下假设：流体是连续、不可压缩的牛顿流体；温度作为被动标量处理，流体的热物性参数均为常数，即：忽略温度对流场的影响。

1.2.1 控制方程

基于上述假设，CIJR内部的流动及混合过程满足以下方程：

连续性方程

$ \nabla \cdot \boldsymbol{U}=0. $

(1)

动量守恒方程

$ \frac{\partial \boldsymbol{U}}{\partial t}+\boldsymbol{U} \cdot \nabla \boldsymbol{U}=-\nabla p+\nu \nabla^2 \boldsymbol{U}. $

(2)

温度输运方程

$ \frac{\partial T}{\partial t}+\boldsymbol{U} \cdot \nabla T=D \nabla^2 T. $

(3)

其中：U为速度矢量；T为温度；P为压强。

1.2.2 SST-IDDES湍流模型

本文采用的SST-IDDES模型基于SST k-ω模型进行改造^[34]，通过引入湍流长度尺度l_IDDES对模型中的湍动能耗散项进行修正，l_IDDES定义为

$ l_{\mathrm{IDDES}}=\widetilde{f_{\mathrm{d}}}\left(1+f_{\mathrm{e}}\right) l_{\mathrm{RANS}}+\left(1-\widetilde{f_{\mathrm{d}}}\right) l_{\mathrm{LES}} . $

(4)

其中：l_LES=C_DESΔ_IDDES和$l_{\mathrm{RANS}}=\frac{\sqrt{k}}{C_\mu \omega}$分别为LES滤波尺度和RANS长度尺度；C_DES=0.78F₁+0.61(1- F₁)和C_μ=0.09是模型系数，F₁为SSTk-ω模型的混合函数；Δ_IDDES是亚格子尺度，定义为

$ \Delta_{\text {IDDES }}=\min \left\{\max \left[0.15 d_{\mathrm{w}}, 0.15 h_{\max }, l_{\mathrm{w}}\right], h_{\max }\right\} . $

(5)

该定义有助于解决对数层不匹配(log-layer mismatch，LLM)问题。其中，d_w是计算点到壁面的距离，h_max=max(Δx, Δy, Δz)是当地网格的最大长度，l_w是垂直壁面方向的网格长度。

式(4)中$\widetilde{f_{\mathrm{d}}}$为混合函数，按下式计算：

$ \widetilde{f_{\mathrm{d}}} =\max \left[\left(1-f_{\mathrm{dt}}\right), f_{\mathrm{B}}\right], $

(6)

$ f_{\mathrm{dt}} =1-\tanh \left[\left(20 r_{\mathrm{dt}}\right)^3\right] . $

(7)

其中，经验混合函数f_B用来实现在壁面距离0.5h_max < d_w < h_max范围内由RANS模型(f_B=1)到LES模型(f_B=0)的快速转换，定义为

$ f_{\mathrm{B}}=\min \left[2 \exp \left(-9 \alpha^2\right), 1.0\right], $

(8)

$ \alpha=0.25-d_{\mathrm{w}} / h_{\max } . $

(9)

式(4)中f_e为经验函数，用于抑制RANS/LES交界区域雷诺应力的过度减小，也有助于解决LLM问题，定义为

$ f_{\mathrm{e}}=\max \left[\left(f_{\mathrm{e} 1}-1\right), 0\right] f_{\mathrm{e} 2}, \\ $

(10)

$ f_{\mathrm{e} 1}=\left\{\begin{array}{cc} 2 \exp \left(-11.09 \alpha^2\right), & \alpha \geqslant 0, \\ 2 \exp \left(-9 \alpha^2\right), & \alpha <0, \end{array}\right. $

(11)

$ f_{\mathrm{e} 2}=1-\max \left(f_{\mathrm{t}}, f_1\right), $

(12)

$ f_{\mathrm{t}}=\tanh \left[\left(C_{\mathrm{t}}^2 r_{\mathrm{dt}}\right)^3\right] \text {, } $

(13)

$ f_{\mathrm{l}}=\tanh \left[\left(C_1^2 r_{\mathrm{dl}}\right)^{10}\right], $

(14)

$ r_{\mathrm{dt}}=\frac{\nu_t}{\max \left(\sqrt{\left(\partial u_i / \partial x_j\right)\left(\partial u_i / \partial x_j\right)}, 10^{-10}\right) \cdot \kappa^2 d_{\mathrm{w}}^2}, $

(15)

$ r_{\mathrm{dl}}=\frac{\nu}{\max \left(\sqrt{\left(\partial u_i / \partial x_j\right)\left(\partial u_i / \partial x_j\right)}, 10^{-10}\right) \cdot \kappa^2 d_{\mathrm{w}}^2} . $

(16)

其中：参数r_dl和r_dt分别代表层流和湍流成分；C_t和C_l是模型常数，取决于所采用的RANS模型。

在实际计算中，根据来流条件的不同，IDDES方法可以分为DDES和WMLES 2个分支。当来流不包含湍流成分时，f_e=0，IDDES采用DDES分支计算，IDDES长度尺度简化为

$ l_{\mathrm{IDDES}}=\widetilde{f_{\mathrm{d}}} l_{\mathrm{RANS}}+\left(1-\widetilde{f_{\mathrm{d}}}\right) l_{\mathrm{LES}} . $

(17)

相反，当来流包含湍流成分时，$\widetilde{f_{\mathrm{d}}}=f_{\mathrm{B}}$，IDDES采用WMLES分支计算，IDDES长度尺度简化为

$ l_{\mathrm{IDDES}}=f_{\mathrm{B}}\left(1+f_{\mathrm{e}}\right) l_{\mathrm{RANS}}+\left(1-f_{\mathrm{B}}\right) l_{\mathrm{LES}} . $

(18)

1.3 网格划分、初始条件和边界条件 1.3.1 网格划分

本文采用ICEM CFD18.0软件对计算域进行结构化网格划分。考虑到撞击区和SRZ附近的流动比较复杂，为了保证计算精度，在划分网格时对2个区域进行了网格加密处理。为了满足不同区域对分辨率的要求，同时将计算成本控制在合理范围，本文设计了3套不同的网格划分方案，如表 1所示。

Δ_i⁺为无量纲网格尺寸，定义为

$ \Delta_i^{+}=\Delta_i \frac{u_\tau}{\nu}(i=x, y, z, ). $

(19)

式中：Δ_i(i=x, y, z)分别代表流向、横向和展向的实际网格尺寸；u_τ为摩擦速度。Δ_i,bulk⁺代表腔室中心网格的无量纲尺寸，Δ_i,wall⁺为壁面边界层内第一层网格的无量纲尺寸。网格方案Ⅰ中的网格数量约为360万，上游重要区域(-0.5 < X^* < 2.0；-1 < Y^* < 1；-1 < Z^* < 1)内的网格比较精密，Δ_i,bulk⁺设置为50，沿流向网格尺寸逐渐增大。网格Ⅱ和网格Ⅲ在网格Ⅰ的基础上对撞击区和SRZ的网格进行局部加密。3套方案对壁面附近的网格适当加密保证Δ_i,wall⁺ < 1，以满足IDDES模型计算时对网格分辨率的要求^[34]。本文中，上标*用于识别无量纲量，3个坐标变量分别使用H_C、L_S、W_C/2进行无量纲化；3个速度分量均用体积速度U_b进行无量纲化；流场温度使用射流温差ΔT进行无量纲化。

1.3.2 初始条件和边界条件

本文计算中主要涉及的边界条件有以下3种：

1) 入口边界2个喷嘴入口设置为速度入口边界条件。入口射流体积速度为U_b=13.49 m/s，对应的雷诺数$R e=\frac{U_{\mathrm{b}} D_{\mathrm{H}}}{\nu}=41\; 000$，其中D_H为水力直径；ν为空气的运动黏度。在均匀来流下，通过直接给定湍流强度的方式定义入口湍流。2个喷嘴入口处射流的温度剖面设置为均匀的，射流温差为ΔT=T₂-T₁=1 K(T₂>T₁)，其中T₁和T₂分别是射流1和射流2的温度。

2) 出口边界反应器的腔室沿流向足够长，在出口处流动达到充分发展状态。

3) 壁面边界喷嘴和腔室的固体壁面均设置为无滑移和绝热边界条件。

初始时刻，反应器内部流体处于静止状态，温度为$T_0=\frac{T_1+T_2}{2}$。

1.4 数值求解方案

本文数值计算基于开源流体计算软件OpenFOAM的平台，在非反应流不可压缩求解器pisoFoam中添加温度输运方程并重新编译，对CIJR内的撞击混合过程进行IDDES模拟。所有方程中的时间导数项均采用二阶隐式Crank-Nicolson格式离散。动量方程中的对流项采用线性-迎风稳定输运格式进行离散，该格式混合中心差分格式和迎风格式，以消除复杂流动中无界中心差分格式引起的非物理振荡。输运方程中的对流项采用二阶Van Leer格式进行离散。动量方程中的扩散项采用二阶中心差分格式进行离散。为了解决压力-速度耦合问题，使用PISO(pressure-implicit with splitting of operators)算法对离散后得到的方程组进行求解。计算过程中，选择的时间步长满足最大局部Courant数小于1。

2 模型验证和网格无关性验证

为验证数值模型的可靠性和网格的无关性，本节将CIJR流场的IDDES模拟结果与Liou等^[20]的实验数据进行详细比较。

图 2(a)和2(b)为采用3种网格方案计算得到的展向中心平面(Z^*=0)上不同流向位置的无量纲湍动能和雷诺应力沿横向Y的分布。其中，网格方案Ⅰ中的网格数量约为360万；网格方案Ⅱ和Ⅲ对撞击区(0 < X^*≤0.5)和SRZ(0.5 < X^* < 1.67)的网格进行了加密，网格数量分别为700万和1 030万。从图中可以看出，在圆顶区(X^* < 0)和远下游区域(X^*≥4.0)，采用3种网格的计算结果均与实验结果吻合良好。在撞击区和SRZ，采用网格Ⅰ计算得到的结果与实验结果之间存在较大偏差；而采用网格Ⅱ和网格Ⅲ计算得到的结果与实验结果接近。当1.1 ≤ X^* < 4.0时，本文模拟结果与实验结果之间的误差相对较大。本文模拟得到的湍动能和雷诺应力剖面均关于Y^* = 0对称，而实验数据不具有这种对称性。作者认为这种差异主要源于样本量的差异。本文及实验结果^[20]均表明：流动在统计上呈现周期振荡特性，主导频率为250 Hz。本文的统计平均量是基于30个振荡周期的流场信息，样本量较大。对于实验数据，Liou等^[20]没有提供相关信息。综上所述，采用网格Ⅱ和网格Ⅲ得到的IDDES结果可以较好地捕捉流场中的大尺度结构特征。为了保证计算结果的准确性并节约计算成本，后续对流场的分析与讨论将采用网格Ⅱ的计算结果。

3 相平均方法

相平均方法是基于Matsumura和Antonia^[35]提出的相平均原理，对相位的定义进行了改进。首先，使用快速傅里叶变换对(X^*, Y^*, Z^*)=(0.5, 0, -0.5)处的横向速度V的时间序列进行频谱分析，发现在频率f=250 Hz处出现明显的峰值，表明流场中存在周期性振荡现象。然后，选择二阶低通滤波器Butterworth对原始信号V进行滤波得到滤波信号V_f，滤波频率设置为250 Hz。根据滤波后的信号V_f定义相位

$ \varphi=\pi\left(\frac{t-t_{1, i}}{t_{2, i}-t_{1, i}}\right), \quad\left(t_{1, i} \leqslant t \leqslant t_{2, i}\right), $

(20)

$ \varphi=\pi\left(\frac{t-t_{2, i}}{t_{1, i+1}-t_{2, i}}\right)+\pi, \quad\left(t_{2, i} \leqslant t <t_{1, i+1}\right) . $

(21)

其中，t_1,i和t_2,i分别对应于V_f=0, dV_f/dt < 0和V_f=0, dV_f/dt>0的时刻，因此区间t_2,i~t_1,i代表流场振荡周期τ的一半(τ=1/f)。以π/2为间隔在相位φ=0~2π之间定义4个相位，并对每个相位进行统计平均。需要注意的是，相平均过程是在原始信号上进行的。

瞬时量ϕ可以分解为拟序相平均分量〈ϕ〉和随机脉动分量ϕ′，拟序相平均分量〈ϕ〉可以进一步被分解为时间平均分量ϕ和拟序脉动分量$\widetilde{\phi}$^[36]，即

$ \phi=\langle\phi\rangle+\phi^{\prime}=\bar{\phi}+\widetilde{\phi}+\phi^{\prime} . $

(22)

计算30个振荡周期进行相位平均，得到每个相位的拟序相平均分量〈ϕ〉_φ为

$ \langle\phi\rangle_{\varphi}=\frac{1}{N} \sum\limits_{i=1}^N \phi_{\varphi, i}. $

(23)

其中：φ代表相位编号，N代表周期个数。

4 结果分析与讨论 4.1 不同区域的流动特性和涡结构特征 4.1.1 自由射流区和圆顶区的流动特性和涡结构特征

图 3为一个周期内不同时刻近场区域的展向中心平面Z^*=0上的无量纲瞬时展向涡量场和二维流线图，其中，正涡量使用红色进行着色，而负涡量使用蓝色进行着色，4个瞬时时刻分别与4个相位φ=0、φ=π/2、φ=π和φ=3π/2对应。值得注意的是，这些流线不是真实的流线，而是绘制在展向截面上的二维流线。从流线图可以看出，圆顶区产生一对反向旋转的涡旋，与Liou等^[20]的实验结果吻合。在腔室中心附近，2股射流发生撞击，射流在撞击产生的压力和涡旋的作用下逐渐弯曲并向腔室下游流动。另外，对展向涡量场的观察发现在2股射流的上游边缘存在2个涡量集中的区域，表明射流剪切层中K-H不稳定性发展导致展向涡的卷起，2个剪切层内的展向涡在强度上基本相当，旋转方向相反。进一步地，从展向涡量场可以观察到2个射流剪切层在撞击面附近相互错开，并且周期性地转换方向，即出现周期性偏转振荡现象。这种现象属于自持振荡，由射流的自身不稳定性和2股射流的撞击不稳定性引起，借助于撞击区压力的周期性变化得以维持，与圆顶区对流动的限制作用密切相关。此外，射流剪切层的自持振荡也伴随展向涡的周期性运动。

4.1.2 撞击区的流动特性和涡结构特征

从图 3中还观察到撞击面围绕腔室轴线(Y^*=0)振荡，即撞击驻点出现周期性偏移现象。φ=0时，撞击驻点位于轴线的下方；φ=π/2时，撞击驻点回到轴线附近；φ=π时，撞击驻点向上方偏移；当φ=3π/2时，撞击驻点再次回到轴线附近。这是撞击流流场中的一个重要现象，国内学者李伟锋及其团队对这种现象进行了大量的实验和数值模拟研究^[37-39]。但本文研究发现的撞击驻点的不稳定性和李伟锋等^[37]对两喷嘴对置撞击流流场进行实验研究时发现的驻点偏移规律不完全相同：他们的实验结果表明当2股射流速度不相等时，动量也不相等，在轴向上迎面发生碰撞时，撞击面垂直于轴线向动量较小的一侧偏移；而本文发现撞击驻点出现周期性偏移现象的原因有2个：一方面，2股射流的一部分会交替地偏转进圆顶区，导致2股射流的动量不相等；另一方面，射流剪切层内的展向涡随着剪切层的偏转振荡周期性地运动，因此撞击驻点的偏移轨迹不在射流轴线上。需要指出的是，李伟锋等^[37]的研究对象是自由撞击流，射流的发展不受固体壁面的限制。

撞击区的另一个重要流动特征是产生了大尺度的流向涡^[31-32]。图 4(a)和4(b)分别展示了由无量纲流向涡量着色的瞬时Q=25等值面和无量纲流向涡量ω_x^*=12等值面。同样地，正涡量和负涡量分别使用红色和蓝色进行着色。由图 4(a)可知，流向涡位于腔室中心附近，以反向旋转的涡对形式存在并沿展向交替排列。

根据前人的研究^[32]，在向下游发展的过程中，同向旋转的流向涡相互靠近并发生合并；在SRZ附近，流向涡发生破裂，在图 4(b)中观察到了这些现象。与Liou等^[20]的实验结果显示的4~5对流向涡不同，本研究识别出来的流向涡对的数量较多。我们认为导致这种差异的原因是Liou等^[20]通过测量撞击区沿(X^*=0.28，Y^*=0)这条线上的一维时均横向速度确定涡结构的个数，而由图 4(a)可知，流向涡在横向上分布不规则，有一些流向涡没有穿过(X^*=0.28，Y^*=0)这条线。

4.1.3 下游发展区的流动特性和涡结构特征

下游发展区的流动特性和涡结构特征尚未得到揭示。为识别下游发展区的拟序涡结构，本文定义了一个相对速度等于U_b的坐标系，并在该坐标系中考察展向中心平面Z^*=0上的相平均流线图，如图 5所示。研究发现，除了SRZ，下游发展区还存在大尺度展向涡，SRZ的2个再附着长度的周期性变化对展向涡的形成有重要影响。以相位φ=3π/2为例，上方剪切层的再附着长度相较于下方剪切层的再附着长度更大，根据连续性条件，腔室轴线下方的流体会先经历强减速，因此轴线两侧的流体之间存在较大的速度差，在两侧流体之间的剪切混合层内，K-H不稳定性逐渐发展并产生一个顺时针旋转的展向涡；之后，上方剪切层再附着于壁面，引起轴线上方的流体减速，从而诱发一个逆时针旋转的展向涡的卷起。在经典K-H不稳定性中，2股流体界面附近的速度导数符号通常是固定的，流动失稳导致同向旋转的涡结构。在本研究中，由于流动在统计上呈现周期振荡特性，剪切层附近的速度导数沿流向呈正、负交替变化，因此，涡结构的旋转方向是顺时针和逆时针交替变化的。在下游发展区，共有4~5对旋转方向相反并沿流向交替排列的展向涡。偏转进下游发展区的射流与腔室轴线不平行，因此2股流体之间的湍流混合层会沿横向方向上下振荡。随着流动的发展，2股射流之间的速度差减小，导致展向涡的横向尺度减小，展向涡的涡核沿流向逐渐偏离腔室轴线，在腔室出口处几乎观察不到展向涡的存在。

4.2 不同区域的混合特性 4.2.1 撞击区流向涡对混合特性的影响

为了研究热量输运和涡旋之间的关系，图 6给出撞击区截面X^*=0.28上的无量纲瞬时温度场和速度矢量场。由图可知，瞬时温度场在腔室中心附近存在明显的波动，表明2股射流之间发生了显著混合。一方面，撞击区的流向涡卷吸环境流体，涡旋边界处的流体渗入涡核，实现混合；另一方面，由图 4(b)可知，随着流向涡向下游发展，旋转方向相同的涡相互靠近并发生合并，涡的空间尺度逐渐增大，因此流向涡卷吸的环境流体经历涡量波动引起的介观混合。但由于流向涡的横向和展向尺度远小于腔室尺寸，流向涡只能改善腔室中心附近的流体混合。

4.2.2 下游发展区展向涡对混合特性的影响

图 7给出展向中心平面Z^*=0上的无量纲相平均温度场。由图可知，在2股流体的剪切混合层中无量纲温度接近于0，表明冷、热流体被展向涡交替地输送到腔室中心并发生混合。此外，由于大尺度展向涡的卷吸拉伸，2股流体之间的接触面积增大，从而促进混合。与2股流体之间的湍流混合层相似，随着展向涡向下游发展，热混合层在横向方向上发生大尺度振荡，冷、热流体的温差迅速减小。因此，大尺度展向涡有效促进了热量输运，显著强化了湍流混合。

4.2.3 相平均热通量分布

图 8(a)和8(b)分别展示一个周期内4个相位沿腔室轴线Y^*=0的无量纲相平均流向热通量分布和横向热通量分布。2个剖面的振荡反映了流场的不稳定性。由图可知，4个相位的横向热通量分量的绝对值峰值位于0 < X^* < 1，在此区域内流向热通量分量较小，表明在撞击区存在强烈的横向热量输运。流向热通量分量的绝对值峰值位于2 < X^* < 3，对应第1个展向涡卷起的位置。随着流向距离的增加，观察到2个热通量分量极值位置的周期性移动，2个热通量分量的振荡范围逐渐变窄，意味着随着展向涡横向尺度的减小，展向涡输运的热量逐渐减少。

5 结论

CIJR流场中存在各种类型的拟序涡结构，这些涡结构的动力学行为对流场的三维流动特性和混合特性有重要影响。本文采用IDDES方法对CIJR内流动特性和混合特性进行研究。在此基础上使用相平均方法将拟序结构从湍流流场中分解出来，通过相平均速度场和瞬时速度场较为全面地分析自由射流区、圆顶区、撞击区和下游发展区的流动特性，揭示了涡结构的生成机制和演化特征，并通过相平均温度场和热通量场分析了涡结构对混合的影响。本文得到的主要结论如下：

1) 在射流出口附近的上游剪切层内，K-H不稳定性的发展诱发了展向涡的产生，2股射流剪切层内的展向涡在强度上基本相当，旋转方向相反。由于射流自身不稳定性和2股射流的撞击不稳定性，射流剪切层发生周期性偏转振荡，这种自持振荡与速度和压力的周期性变化以及PRZ涡流的周期性运动密切相关。此外，射流偏转导致2股撞击射流的动量不相等，剪切层的自持振荡伴随展向涡的周期性运动，进而导致撞击驻点向动量较小的一侧偏移，且偏移轨迹不在射流轴线上。

2) 在下游发展区，由于SRZ的2个再附着长度不相等，腔室轴线两侧的流体存在速度差，导致K-H不稳定性在两侧流体之间的剪切层内发展并卷起大尺度的展向涡。相邻展向涡的旋转方向相反并沿流向交替排列，随着流动向下游发展，展向涡的横向尺度减小。

3) 流向涡在撞击区引起了强烈的横向热量输运，由于流向涡的展向和横向尺度较小，热量输运主要发生在腔室横向中心附近。下游发展区的大尺度展向涡对环境流体的卷吸导致2股流体之间的接触面积增大, 有效促进了热量输运，显著强化了湍流混合。因此，探讨撞击流反应器内涡结构的生成和演化机制对于优化几何构型和运行环境实现强化混合的工程应用目标具有重要的理论指导意义。

4) 模拟结果和实验结果能够很好地吻合，说明IDDES方法能够可靠有效地模拟CIJR内部流动。