VORTEX STRUCTURE ANALYSIS OF VORTEX RING COLLISION PROCESS BASED ON DIRECT NUMERICAL SIMULATION
-
摘要: 涡环间的碰撞涉及复杂的多层级涡结构拓扑关系和涡量转移机理. 深入研究该模型, 有助于揭示涡结构之间相互作用过程所包含的机理. 区别于现有研究一般专注于涡环小尺度次级结构的瞬态演化, 本研究使用基于一种改进的SIMPLE算法的直接数值模拟, 结合整体动态演变与其涡结构细节对涡环碰撞模型进行分析. 基于高分辨率网格, 分析了在不同雷诺数下涡环碰撞后由涡重联分裂出的次级涡结构. 研究表明, 涡环碰撞过程伴随着复杂的涡结构拓扑变化, 如高涡量区域的偏转及涡丝间的此消彼长. 在碰撞时, 环本身的不稳定性导致了局部接触后分裂出沿方位角规律分布的次级涡环. 而次级结构本身随流场雷诺数的增加呈现出从无到有、从近似圆环状到湍流态的转变. 最后, 对多雷诺数下的涡拟能曲线进行横向比较以反映碰撞对流场的整体涡流强度影响. 结果表明, 雷诺数的提高对涡拟能的调控整体表现为初始时刻耗散速率降低; 因碰撞而抬升的峰值提高; 峰值对应的时间点先延后再前移. 本研究基于直接数值模拟展示了涡结构演化的偏转、重联、湍流化等特征, 揭示了雷诺数对涡环碰撞的部分调控机制, 对复杂流动的涡系演化机理研究有启示作用.Abstract: Vortex ring collision involves intricate multi-level topological relationships of vortex structures and mechanisms of enstrophy transfer. A thorough investigation of this model contributes to the elucidation of the mechanisms involved in the interaction processes among vortex structures. Distinguishing itself from existing studies that typically focus on the transient evolution of small-scale secondary structures within vortex rings, this research employs direct numerical simulation utilizing an enhanced SIMPLE algorithm. It integrates the overall dynamic evolution with detailed analysis of vortex structure to investigate vortex ring collision models. Utilizing a high-resolution grid, we analyzed the secondary vortex structures spawned after vortex ring collisions at different Reynolds numbers. The study reveals that the vortex ring collision process is accompanied by complex topological changes in vortex structures, such as the deflection of high enstrophy regions and the elongation between vortex filaments. During collisions, the instability of the rings themselves leads to the formation of secondary vortex rings distributed according to azimuthal angles after local contacts. The secondary structures undergo a transition from non-existence to existence and from an approximately circular ring shape to a turbulent state with increasing Reynolds numbers in the flow field. Concerning the enstrophy profiles reflecting the overall flow field, an increase in Reynolds number generally results in a reduction of the initial moment of dissipation rate, an elevation in the peak associated with collision-induced uplift, and a shifting forward of the peak time corresponding to the elevation. Overall, variations in Reynolds number significantly influence the vortex flow intensity and evolution of vortex structures at various stages of vortex ring collision. This study demonstrates the deflection, reconnection, and turbulization features of vortex structure evolution based on direct numerical simulations. The partial regulation mechanism of Reynolds number on vortex ring collision is revealed, which is instructive for the study of the evolution mechanism of vortex systems in complex flows.
-
Keywords:
- vortex ring collision /
- vortex reconnection /
- instability /
- enstrophy /
- direct numerical simulation
-
引 言
涡环是流体形成封闭环状涡管的现象, 广泛存在于生物运动[1]、飞机尾流[2]和工业切割[3-4]等多种天然和人工环境. 多个涡环在流场中交汇引发的涡环碰撞现象在近年来由于其复杂的能量转移和流动结构演化而引起了海内外学者的关注[5].
在实验研究方面, 涡环及其碰撞实验曾揭示许多独特的机制[6]. Lim等[7]与Chu等[8]于20世纪90年代先后完成了碰撞试验: 采用染色方法比较细致地观察到了次级涡结构, 并对碰撞过程中短波不稳定性及不同雷诺数下的湍流云进行了讨论. Mckeown等[9]运用PIV与数值模拟对涡环碰撞中的级联问题进行深入探究, 首先提出了涡环碰撞的结构发展过程: 涡环首先会在局部出现帐篷型的扭曲, 变为涡旋薄片, 然后成为次级涡结构, 进而演变成3级的涡结构, 并最终发展至完全耗散.
相对于实验研究, 数值仿真可以更精细地捕捉流场的结构与动能、动量输运过程. Cheng等[10]采用格子玻尔兹曼方法进行了直接数值模拟. 结果表明, 在特定的雷诺数下, 具有不同环径比的两个涡环会发生相互滑移, 随后继续前进并逐渐分离. Guan等[11]采用可压缩的有限体积方法进行大涡模拟. 结果表明, 直接碰撞与涡环本身对流体的夹带作用共同主导着环的核心涡度变化.
在多数涡环碰撞的数值仿真中, 为确保计算结果与实验观测相匹配, 通常需要引入人工扰动[12]. 目前研究中采用的扰动模型大多源于对Crow不稳定性[13]和椭圆不稳定性[14-15]的研究. Crow不稳定性具体表现为涡重联现象: 两条细长的涡线展现出对称且类似于正弦波的长波不稳定性, 这导致涡线在诱导效应下互相交互, 间隔性地形成一系列小涡环, 如图1中Van Hooft等[16]所得计算结果. 椭圆不稳定性是由于流场中的扭曲应变场使涡线从圆形截面变为椭圆形, 从而引发的短波不稳定性[17].
![]()
不稳定性理论和人工扰动对于涡环碰撞的数值模拟有着重要影响. Cheng等[18]引入扰动模型后观察到, 涡核较厚的涡环在中等雷诺数下会出现方位角扰动的重排列; Mishra等[19]通过采用级数形式的正弦波扰动, 并结合速度模态分析等技术进行研究, 其涡量范数的计算结果如图2所示. 研究发现, 在不同的雷诺数下, 两种不稳定性的效果并不相同. 当涡环的环径比达到一定水平时, 其涡旋行为不能仅被归纳为类Crow或类椭圆的不稳定性, 还需要考虑涡的结构变化.
![]()
基于不稳定性理论的数值模拟已取得诸多进展, 但计算能力的局限使得参数分析和机理研究还面临挑战. 高分辨率模拟为小尺度涡结构分析提供了丰富的速度场和涡量场数据, 对于一些涡结构交互细节如奇点问题[20]、级联问题[21]乃至带有尾迹的“剥落”等现象[22-25]的自旋涡环机制能够提供重要参考. 这些模型对小结构的分析精度要求推动了研究者追求更高的网格分辨率和仿真精度. 高分辨率直接数值模拟可更精确地捕捉涡环碰撞的细节与流场结构, 揭示涡环碰撞的深层物理机制[26].
在小尺度的涡结构细节分析上, 直接数值模拟(DNS)可更精确地捕捉涡环碰撞的细节与流场结构, 揭示涡环碰撞的深层物理机制. Kollmann等[26]对柱坐标下的碰撞模型进行仿真, 通过涡结构形态、能谱等方法对不稳定的发展进行研究. Hussain等[27]利用基于伪谱算法的直接数值模拟, 对高雷诺数下碰撞产生的二级及以上的涡结构重联演变进行分析, 并对其涡奇点问题进行理论模型构建.
有别于上述研究一般更注重于研究涡结构在瞬态下的发展, 缺少整体的演变发展, 本文将针对涡环碰撞的整体演变进行研究, 主要关注动态演化结合高涡量区的小尺度结构行为.
本文基于一种改进的SIMPLE算法(semi-implicit method for pressure-linked equations)的直接数值模拟对涡环碰撞的过程及其演化进行研究, 特别是因参数变化导致的涡结构的变迁. 在第1部分, 本文将详细描述数值方法, 包括控制方程选取与求解方法、涡环的初始速度场的构建公式以及数值验证. 第2部分将呈现涡环碰撞的数值结果和从不同视角的可视化效果, 其中第1小节将探讨低雷诺数下涡环碰撞的涡结构变化; 第2小节重点分析中等雷诺数下涡结构的演化, 深入探讨三维涡系结构拓扑结构的变化、碰撞过程中的涡量等值线的转变, 并详细描述涡环伸展阶段的涡结构行为; 第3小节将探讨高雷诺数下次级涡环如何从清晰的结构逐渐转变为湍流态; 第4小节将基于涡拟能的趋势, 总结雷诺数提高过程中的标志性特点. 第3部分将对全文进行总结.
1. 数值方法
1.1 流场构造
本研究对两个初始时刻平行、同轴、速度环量相同且异号放置的涡环进行仿真. 求解采用直接数值模拟以三维不可压缩Naiver-Stokes方程作为控制方程, 即
$$ \left. \begin{split} & \nabla \cdot {\boldsymbol{u}} = 0 \\ & \rho \frac{{\partial {\boldsymbol{u}}}}{{\partial t}} + \rho \left( {{\boldsymbol{u}} \cdot \nabla } \right){\boldsymbol{u}} = - \nabla p + \mu {\nabla ^2}{\boldsymbol{u}} \end{split} \right\} $$ (1) 其中, $ t $为时间, ${\boldsymbol{ u}} $为速度, $ p $为压力, $ \rho $为密度, $ \mu $为动力黏度.
流场计算区域设置为笛卡尔坐标系内的长方体区域
$$ \left[ { - \frac{{{L_x}}}{2},\frac{{{L_x}}}{2}} \right] \times \left[ { - \frac{{{L_y}}}{2},\frac{{{L_y}}}{2}} \right] \times \left[ { - \frac{{{L_z}}}{2},\frac{{{L_z}}}{2}} \right] $$ 采用均匀六面体网格进行离散, 边界均采用固壁边界.
在无扰动情况下, 涡环碰撞的初始状态由初始速度环量(涡量)分布确定, 如图3所示. 其中, $ {R_0} $为涡环初始半径, $ {S_0} $为涡环初始间距.
![]()
图 3 无扰动时涡环初始状态的特征量Figure 3. Characteristic quantities of the initial state of the vortex ring in the absence of perturbations单个放置于初始涡量场分布为
$$ {\boldsymbol{\omega}} \left( \sigma \right) = \frac{{{\varGamma _0}}}{{\text{π} \sigma _0^2}}{{\mathrm{e}}^{ - {{\left( {\frac{\sigma }{{{\sigma _0}}}} \right)}^2}}}{{\boldsymbol{e}}_n} $$ (2) 其诱导速度场分布为
$$ {{\boldsymbol{u}}_0}(\sigma ) = \frac{{{\varGamma _0}}}{{2\text{π} \sigma }}[1 - {{\mathrm{e}}^{ - {{\left( {\frac{\sigma }{{{\sigma _0}}}} \right)}^2}}}]{{\boldsymbol{e}}_T} $$ (3) 其中, $ {\sigma _0} $为涡环的特征内径, $ {\varGamma _0} $为涡环的特征(速度环量)强度, $ {{\boldsymbol{e}}_n} $为该点沿涡环中线的单位矢量, $ {{\boldsymbol{e}}_T} $为该点围绕环心的速度环量的单位切向量. 对于双涡环模型的初始场, 则是将两个涡环各自的诱导速度场叠加.
相关文献表明[28-30], 初始涡环变形引发的扰动对碰撞过程中次级涡环生成有着直接影响. 本研究中, 涡环变形由以下公式表示
$$ R = {R_0}\left[ {1 + \varepsilon \cos \left( {n\varphi + {\varphi _0}} \right)} \right] $$ (4) 其中, $ \varepsilon $为扰动强度, $ n $为扰动波数, $ {\varphi _0} $为扰动相位.
引入以下无量纲变量与参数
$$ T = \frac{{{\varGamma _0}t}}{{{R_0}^2}},{Re} = \frac{{\rho {\varGamma _0}}}{\mu },\eta = \frac{{{\sigma _0}}}{{{R_0}}},{\boldsymbol{\varOmega }} = \frac{{{\boldsymbol{\omega }}{R_0}^2}}{{{\varGamma _0}}} $$ (5) 其中, $\eta $为环径比, 用于衡量涡核相对粗细程度. 无量纲变量引入涡环半径${R_0}$与速度环量${\varGamma _0}$以衡量涡环参数对模型演变的影响.
1.2 计算方法与数值验证
本文基于改进的有限差分法与SIMPLE算法求解三维不可压缩Naiver-Stokes方程. SIMPLE算法优势包括数值稳定性高、适用范围广泛等. 然而, 该算法也存在如可能引入数值扩散、迭代步骤多导致计算成本高等问题. 对于计算效率, 本文对SIMPLE算法中的压力修正过程进行改进, 将传统的迭代收敛方法改为基于快速傅里叶变换对压力泊松方程矩阵进行单次求解, 提高计算效率.
改进的算法在高精度的流场仿真基础上, 相较同类型算法的计算效率有显著提高. 改进算法的对流项采用守恒型差分, 扩散项采用中心差分离散. 为实现对涡环问题的直接数值模拟, 本研究采用了一种改进的基于FFT算法的Poisson方程快速求解方法. 该方法将基于中心差分格式离散的压力Poisson方程转换为多个独立子方程以实现并行计算, 并由FFT算法降低变换过程的计算复杂度. 具体而言, 该算法分为以下5步进行.
(1) 将基于下式中心差分的压力Poisson方程
$$\begin{split} & \frac{1}{{\Delta {x^2}}}\left( {{p_{i + 1,j,k}} - 2{p_{i,j,k}} + {p_{i - 1,j,k}}} \right) + \\ &\qquad \frac{1}{{\Delta {y^2}}}\left( {{p_{i,j + 1,k}} - 2{p_{i,j,k}} + {p_{i,j - 1,k}}} \right) + \\ &\qquad \frac{1}{{\Delta {z^2}}}\left( {{p_{i,j,k + 1}} - 2{p_{i,j,k}} + {p_{i,j,k - 1}}} \right) = {f_{i,j,k}} \end{split} $$ (6) 对应的线性方程组$ {\boldsymbol{Ap}} = {\boldsymbol{f}} $进行子特征分解, 即$ {\left[ {\boldsymbol{A}} \right]_{i,j}} = {\boldsymbol{Q}}{\left[ {\boldsymbol{\varLambda}} \right]_{i,j}}{{\boldsymbol{Q}}^{\mathrm{T}}} $ ($ {\boldsymbol{Q}} $为块矩矩阵的归一化特征向量), 并计算右端项的变换$ {\left[ {\hat {\boldsymbol{f}}} \right]_{i,j}} = {{\boldsymbol{Q}}^{\mathrm{T}}}{\boldsymbol{f}} $;
(2) 重新排列矩阵与右端向量下标, 空间坐标$ \left( {x,y,z} \right) $顺序重排为$ \left( {z,x,y} \right) $并对重排后矩阵的子特征分解与右端项变换;
(3) 再次重排方程, 使空间坐标顺序重排为$ \left( {z,y,x} \right) $, 再次进行子特征变换后方程转变为一系列彼此独立的3对角子方程;
(4) 基于Thomas等低复杂度方法并行求解所有独立的子方程;
(5) 依次对解进行步骤3、步骤2和步骤1的逆变换, 获得在原空间中空间坐标排列顺序为$ \left( {x,y,z} \right) $的解.
一方面, 上述方法可以将Poisson方程由全局求解分解为独立子方程求解以提高计算局部性; 另一方面, 上述计算过程中的子特征分解及其逆变换由于方程的结构可以由FFT算法加速其中的主要过程. 因此, 对于SIMPLE算法求解过程, 该算法显著降低计算的复杂度.
其中, 下文涉及的动能与涡拟能如下
$$ {E_{{u}}} = \frac{1}{2}\iiint\limits_v {{\boldsymbol{u}} \cdot {\boldsymbol{u}}{\mathrm{d}}v},{E_\varOmega } = \frac{1}{2}\iiint\limits_v {{\boldsymbol{\varOmega }} \cdot {\boldsymbol{\varOmega }}{\mathrm{d}}v} $$ (7) 其中${E_{{u}}}$ 为流场动能, ${E_{{\varOmega}} }$ 为流场涡拟能. 为统一尺度、集中讨论统计量的演化趋势, 采用归一化进行统计量处理. 归一化的动能与涡拟能如下
$$ {E_{{u}}}^* = \frac{{{E_{{u}}}}}{{{E_{u0}}}},{E_{{\varOmega}} }^* = \frac{{{E_{{\varOmega}} }}}{{{E_{\varOmega 0}}}} $$ (8) 其中, ${E_{{{u}}0}}$和${E_{{{\varOmega}} 0}}$分别为初始时刻动能与涡拟能. 后文出现的涡拟能与动能等均为归一化处理后的结果.
在涡环碰撞模型下, 对该算法的独立性进行涡拟能峰值的校验, 如表1所示.
表 1 涡拟能峰值的独立性校验Table 1. Independent verification of the peak value of enstrophydh dt 4.0 × 10−4 2.0 × 10−4 1.0 × 10−4 5.0 × 10−5 1/60 8.364 13.6964 8.6727 9.3375 1/80 13.6964 8.287 7.7339 8.9026 1/100 31.4801 8.5343 7.4551 8.2894 1/128 — 12.0731 6.5894 6.8445 1/150 — 19.575 5.3496 6.1625 由表1计算结果显示, 涡拟能在网格加密和时间步长缩短的趋势下, 整体呈现收敛趋势. 综合考虑计算效率与涡结构细节, 选取1/128的网格密度与1.0 × 10−4网格步长.
本文将仿真结果中涡结构等值面与Lim等[7]基于染料追踪方法的实验结果进行对比. 初始条件设置流场雷诺数${Re} = 1600$, 圆环初始距离${S_0} = 2{R_0}$, 环径比$\eta = 0.25$, 扰动波数幅值$ \varepsilon = 0.2 $, 波数$ n = 18 $, 相位$ {\varphi _0} = 0 $. 等值面如图4所示.
![]()
图 4 涡结构变化结果对比$({{Re}} = 1600)$Figure 4. Comparison between the vortex transformations $({{Re}} = 1600)$图4涡结构变化与文献符合情况良好, 对主环演变姿态以及次级涡环的有着准确的复现.
为基于统计量进一步验证本文所采用数值方法的准确性与数值结果的可靠性, 对雷诺数${Re} = 1000$下的涡环与Cheng等[10]的动能和涡拟能仿真结果对比如图5和图6所示. 计算条件设置圆环初始距离${S_0} = 2{R_0}$, 环径比$\eta = 0.1$. 动能结果如图5所示, 涡拟能的对照结果如图6所示.
![]()
图 5 动能发展结果对比(${Re} = 1000$)Figure 5. Comparison between the kinetic energy developments (Re = 1000)![]()
图 6 涡拟能发展结果对比$({Re} = 1000)$Figure 6. Comparison between the enstrophy developments (Re = 1000)由图5和图6的结果显示, 本算法在与Cheng等[10]有关动能与涡拟能的校验中均有着良好的贴合.
综上, 本文所采用的数值离散与求解方法与文献中数值模拟和实验结果均吻合良好, 且本研究所采用的网格分辨率可实现对涡环碰撞的瞬态涡结构的高精度捕捉.
2. 结果与分析
涡环碰撞随雷诺数的变化而在碰撞过程中呈现不同的姿态. 本节将展示其具体的变化过程, 包括涡系结构的拓扑变化及各种数值的走向.
流场的整体设置以初始涡环半径${R_0}$为单位, 计算域设置为${L_x} \times {L_y} \times {L_z} = 5{R_0} \times 25{R_0} \times 25{R_0}$, 网格量对应为$128 \times 640 \times 640$, 以尽量减少方体计算域边界对涡环碰撞充分发展后的结构的阻碍作用. 环的初始相隔距离${S_0} = 3{R_0}$, 扰动幅值$ \varepsilon = 0.05 $, 扰动波数$ n = 16 $. 同时, 根据Mishra等[19]的研究, 碰撞向湍流态过渡的雷诺数阈值总体与其初始环径比$ \eta $成正相关, 下文固定环径比$ \eta = 0.25 $讨论其涡结构状态.
初始预设的扰动波数$ n $直接影响碰撞的次级结构, 而根据Lim等[7]的实验结果, 其环境下充分发展后的次级涡环数一般为15 ~ 20个. 以下设置$ n $分别为15, 16, 17和18, 对其计算结果进行涡拟能与涡量范数$|{\boldsymbol{\varOmega }}|$分布(下文简称涡量)比较, 如图7所示.
![]()
图 7 不同扰动波数$ n $下涡量分布Figure 7. Vorticity under different wave numbers $ n $ of perturbation根据图7的涡量分布显示, 在涡环发生次级涡环演变前, 其涡量发展姿态基本同步. 而在次级涡环基本形成后, 展现出波数$ n $越小, 发展提前, 反之发展滞后的趋势. 考虑到发展速度差异有限, 总体演变并不因波数差异脱离碰撞的一般性规律, 故选取对称性强, 间隔角度特殊性高的波数$ n = 16 $作为下文分析波数. 同时, 设置扰动幅值$ \varepsilon = 0.05 $.
2.1 低雷诺数下涡环碰撞分析
选取雷诺数${Re} = 300$下$T = 30$时碰撞涡量分布如图8所示. 在此雷诺数下, 充分发展时, 涡环有一定程度的主环扩大与环径比减小. 涡环的不稳定性显现, 因此不出现次级结构, 只表现一定程度的相互接近, 并直至耗散完全.
由于涡环初始位置的对称性, 现将$x = 0$所在平面设为涡环碰撞的中面, 下文同理. 此雷诺数下的碰撞中面的涡量等值线如图9所示. 环内部涡量等值线近圆, 基本不表现不稳定性.
![]()
图 8 $T = 30$时雷诺数${Re} = 300$下涡量分布Figure 8. Vorticity distribution under Reynolds number ${Re} = 300$ at $T = 30$![]()
图 9 时间$T = 30$时碰撞中面涡量等值线Figure 9. Collision midplane vorticity contours at time $T = 30$2.2 中等雷诺数下涡环碰撞分析
雷诺数${Re} = 1600$下碰撞过程的涡量分布如图10所示. 该雷诺数下, 主环与次级结构的涡量分布均表现规则的形态, 可作为解释涡环碰撞演变的典型算例.
在此雷诺数下, 涡量分布呈现典型的阶段性态势: 初期, 由图10(a)所示, 涡环初步接触, 不稳定性不明显, 整体呈现平滑圆环状. 图10(b)和图10(c)时, 随着进一步接触, 环的相互作用增强, 各自的诱导速度使对方主环的半径增大, 内径减小, 主环从光滑单环发展为波浪状细环. 继续发展至$T = 24$, 如图10(d)和图10(e)所示, 双环的波浪型特征相互作用, 开始形成次级结构, 原主环部分进一步细长, 主环开始被次级结构所替代. 当发展到$T = 37.5$, 如图10(f)所示, 此阶段下, 主环基本完成了解体, 涡量转移至次级结构继续发展, 形成形态清晰的小涡环, 涡环的数量与方位角上初始余弦的波数相同.
对涡环伸展阶段演变过程, 下文结合涡量分布及碰撞中面涡量等值线进行阶段性分析.
分析图11涡量等值线得, $T = 8$时, 中面的涡量剖面与环的形状大致一致. 当两个环处于早期接触阶段时, 彼此相互作用尚浅, 使得环的半径尚未发生明显扩大. 在这个阶段, 环径比$ \eta $较大, 环表面的形状相对平滑. 具体表现为: 环表面的涡量等值线沿着环的表面连续分布, 呈现平缓的近似圆形, 表明此阶段环仍然是一个相对连续的结构. 与环的表面不同, 在环接触的中面, 其位于环内外径之间的涡量起伏较为明显.
![]()
图 11 时间$T = 8$时碰撞中面涡量等值线Figure 11. Collision midplane vorticity contours at time $T = 8$图11从几何特征上, 碰撞中面的局部峰值区域涡量大小分布沿着极坐标方向呈现间隔规律的峰状排列, 局部峰值区域对应的涡量等值线各自包围封闭. 现以碰撞中面的涡量极大值点定义为局部峰值区的中心位置, 由距平面中心的半径${r_c}$与极坐标角度$ \alpha $决定, 如图11的局部放大图示.
局部峰值区的数量与扰动的波数相同. 同时, 局部峰值区的中心位置与扰动在极坐标系统下的角度相关, 对应扰动公式(4)中$R = {R_0}$处.
随着两者的进一步接近, 环径比显著减小, 相互作用的强度也随之增强. 如图12(a)所示, $T = 14$下涡环内部的高涡量区域随着环径比的变化而逐渐显现在表面上, 呈现为沿方位角方向的柱状结构. 整体来看, 环的形状开始显露出类波浪状的几何特征: 高涡量区域对应两个主环在几何上相对贴合的部分. 而在两侧的涡结构还未完全接触, 导致环之间形成一定的空隙.
在$T = 14$下的碰撞中面涡量等值线如图12(b)所示, 与前一阶段相比, 不稳定性下的类波浪状几何特征持续翻转, 最终随主环的接近相对稳定下来. 其峰值区中心角度变化前后对比如图12(c)所示. 其中, 为将中心位置对齐以观察角度差异, 将中心半径做如下处理
$$ {r_{c1}} = {r_c} - {{\bar r_c}} + 1.0 $$ (9) 其中, $ {{\bar r_c}} $表示将本时刻下的峰值区半径取平均值. 由中心位置变化结果显示, 环内部的局部峰值区经历了1/2个扰动周期的错位. 也即, 在接近后的短时间内, 高涡量区域发生了位置上的重构. 演变后的局部峰值区的几何位置仍然与扰动的几何特性相关联, 但其位置从原本对应扰动公式中$R = {R_0}$处变为对应扰动公式(4)中$R = (1 \pm \varepsilon ){R_0}$处.
当$T = 16$时, 如图13(a)所示, 随着涡环半径的进一步增大, 环径比进一步减小, 环间空隙增大. 根据涡量等值线, 如图13(b)所示, 处于环面位置的等值线相对前时间步进一步向涡环内部的局部峰值区贴近, 环靠近内径的等值线波动加剧.
![]()
图 13 时间$T = 16$时涡量分布Figure 13. Distribution of vorticity at different Reynolds numbers for $T = 16$如图14(a)所示, 当$T = 24$, 主环继续向外扩张, 环径比进一步减小. 同时在主环的左右两侧出现细涡丝. 细涡丝横跨主环所在平面, 涡量高于空隙区, 被称为“桥”涡丝[21], 将两个主环连接在一起. 这些细涡丝构成了次级环初始的涡量聚集区, 涡丝中心所含有的涡量大小与局部峰值区中心的涡量值相当.
如图14(b)所示, 随着主环半径的扩大, 环内径进一步被压缩, 导致局部峰值区的形状变得扁长, 峰值区已经基本将原本的主环分隔开来. 与此同时, 峰值区两端出现明显的翘起, 对应于三维可视化中的“桥”涡丝.
当$T = 30$时, 碰撞中面涡量等值线如图15所示. 原局部峰值区仍保持着扁长的形状, 但其高峰值区域已在中间部分发生了断裂. 断裂的部分与原峰值区左右相邻的同样断裂的部分相结合.
![]()
图 15 时间$T = 30$时碰撞中面涡量等值线Figure 15. Collision midplane vorticity contours at $T = 30$$T = 30\sim 37.5$的涡量演变如图16所示. 空隙区通过“桥”涡丝的伸长和包围初步形成一个完整的环状结构, 即次级涡环, 并集中了剩余的涡量, 开始向外运动. 次级涡环生成的位置大致对应扰动公式(4)中$R = {R_0}$处.
![]()
图 16 时间$T = 30$ ~ $ 37.5$涡量分布演变Figure 16. Evolution of vortex distribution from time $T = 30$ ~ $37.5$如图16所示, 结合$T = 30$, $33.75$, $ 37.5$的演变, 次级涡环生成的涡结构经历了小结构上的演变过程: 起初, 两条的“桥”涡丝占据半圆包围空隙区, 成为次级涡环区的涡量集中处. 在演变过程中, 新生成两条沿主环方位角的涡丝, “桥”涡丝的涡量值相比新生成的涡丝占比下降1, 涡丝逐渐消失. 新生成的涡丝重新将空隙云团包围, 与原本的 “桥”涡丝一同构成次级环.
最后, 由于次级涡环已不再涉及更多的涡结构相互作用, 涡量逐渐下降至终, 标志着碰撞演变过程的结束.
2.3 高雷诺数下涡环碰撞分析
雷诺数${Re} $取3000, 4000和5000时的$T = 30$时刻涡量分布, 如图17所示. 由图可以看出, 涡环演变中涡结构整体随雷诺数的提高呈现出湍流化的趋势.
![]()
图 17 $T = 30$时不同雷诺数下涡量分布Figure 17. Vorticity distribution at different Reynolds number for $T = 30$$T = 8$不同雷诺数下碰撞中面涡量等值线如图18所示. 在前期相互作用较小, 半径刚开始扩大的自由运动阶段, 随着雷诺数的增加, 环整体结构保持完整. 环内部的高涡量区域呈现出类似中等雷诺数下的周期性分布, 如图18, 与扰动的几何特征具有较高的一致性. 涡环的自诱导运动在较高雷诺数下对涡环内部的扰动影响相对有限, 并不能完全破坏其扰动. 在这3种雷诺数下, 与雷诺数${Re} = 1600$下的等值线如图13不同, 碰撞中面环面上的涡量等值线紧贴局部峰值区, 更为凌乱.
![]()
图 18 T = 8不同雷诺数下碰撞中面涡量等值线Figure 18. Collision midplane vorticity contours at different Reynolds number for T = 8$T = 16$时不同雷诺数下涡量分布如图19所示, 随着半径的扩大, 涡量沿环的分布开始呈现差异性: 雷诺数${Re} = 3000$时, 高涡量区域的特征与中等雷诺数如图13(a)的柱状涡丝不同, 呈现出多段的丝状涡结构. 涡量丝连接着相邻两个低涡量值的空隙区; 雷诺数${Re} = 4000$下, 如图19(b), 相比图19(a)完整且平滑的涡量丝, 涡结构呈现断裂、混乱的短条涡丝; 雷诺数${Re} = 5000$下, 对应的高涡量区域已经破裂, 呈现出散点状的涡量分布.
![]()
图 19 $T = 16$时不同雷诺数下涡量分布Figure 19. Distribution of vorticity at different Reynolds number for $T = 16$$T = 30$时不同雷诺数下的涡量分布如图20所示. 当主环完全发展且次级涡环生成时, 高雷诺数下的次级结构不同于由垂直于主环与沿主环各两段的涡丝构成的典型特征, 如图15所示. 在雷诺数${Re} = 3000$时, 如图20(a)所示, 次级结构依然呈现出较明显环状的特征, 但构成次级环的涡线众多, 相互缠绕、扭曲连接, 使得单个次级涡环呈现出较为浑厚的圆团状. 在雷诺数${Re} = 4000$下, 如图20(b)所示, 次级结构为夹带着散乱的丝状的涡量湍流云, 不为清晰环状. 在此雷诺数下, 湍流云表面分布有分散的高涡量区域, 整体上仍然呈现出明显的团状结构. 雷诺数${Re} = 5000$下, 如图20(c)所示, 次级结构在充分发展后完全为湍流云团, 云团本身形态扭曲, 整体为不规则的团状雾区. 此雷诺数下, 在发展的湍流云团之间, 还夹杂着残留的涡结构.
![]()
图 20 $T = 30$时不同雷诺数下涡量分布Figure 20. Distribution of vorticity at different Reynolds number for $T = 30$高雷诺数下由于没有形成明确的次级涡环, 也没有结构间的相互作用, 次级结构以这种形态耗散完毕.
2.4 不同雷诺数下涡拟能演变特征
涡拟能可以反映涡环碰撞对流场的整体涡流强度的影响. 现选取雷诺数${Re} $分别为300, 1200, 2100, 3000和3900下归一化涡拟能全场积分结果如图21所示. 从图中可以看出, 整体而言, 涡拟能的升降与碰撞阶段联系密切. 在初始时刻, 由于两环相距较远, 自由运动使得环涡量被消耗, 涡拟能下降. 在两环接近, 相互作用增强后, 涡核涡量增强, 涡拟能随之爬升. 最后, 在碰撞后期, 相互作用基本结束后, 耗散作用使得涡拟能下降至终.
通过涡拟能曲线之间的比较可以发现, 雷诺数的变化对涡拟能曲线的影响明显. 接下来通过提取曲线特征信息进一步分析.
碰撞发展过程中涡拟能峰值高度随雷诺数的变化如图22所示. 整体而言, 碰撞导致的涡核强度增大随雷诺数的提高而持续提高, 且在测试范围内, 雷诺数${Re} = 2500$以上曲线增长速率比2000以下更高, 代表高雷诺数作用于碰撞涡核演变强度变化会更明显.
碰撞发展过程中涡拟能峰值高度所对应时间随雷诺数的变化如图23所示. 雷诺数${Re} $ ≤ 2700时, 其峰值对应时间随雷诺数提升而延后, 但之后, 则又会前移.
![]()
图 21 不同雷诺数下涡拟能在时间上的变化Figure 21. Variation of normalized vortex fitting energy in time for different Reynolds numbers![]()
图 22 涡拟能峰值高度随雷诺数的变化Figure 22. The variation of peak enstrophy height with Reynolds number![]()
图 23 涡拟能峰值对应时间随雷诺数的变化Figure 23. The time position of peak enstrophy with Reynolds number3. 结 论
本研究利用基于一种改进的SIMPLE算法的DNS数值模拟, 探究了不同条件下涡环碰撞全过程, 并进行了定量与定性分析. 在参数分析上, 雷诺数作为主要讨论的流场参数, 对碰撞演变有明显的影响. 涡环碰撞模型的参数包括主涡环的半径、内径及速度环量, 以及两个环初始时的相对距离. 为模拟实验中次级结构的形成, 碰撞模型引入余弦形式的扰动, 其参数包括余弦扰动的波数、振幅和相位.
研究聚焦于不同雷诺数下涡环碰撞过程的形态演变, 分为高、中和低3种雷诺数情况. 在低雷诺数下, 涡环碰撞只表现为一定程度的主环扩大. 主环间不出现明显的涡重联现象, 也不出现清晰的次级结构. 在中等雷诺数下, 涡环碰撞能最清晰地呈现涡环碰撞过程中的涡结构细节. 在此情况下, 涡环表面附近的涡量等值面呈现类环状, 连续且相对平滑. 环内部涡量沿方位角方向存在起伏, 与扰动设置密切相关. 当两个环相距较近时, 相互的诱导作用使得两个环的半径扩大, 涡核变得细长. 此时环的不稳定性逐渐显现在表面, 主环的接触区域之间形成空隙区. 随着碰撞发展, 主环的高涡量区域会向两端空隙区转移, 最终从中间断裂. 空隙区将生成细小的涡丝, 并逐渐形成环形, 次级环生成. 最后, 主环解体, 次级环径向往外运动, 并最终耗散. 在高雷诺数情况下, 在主环碰撞扩大之前, 作用于主环的扰动并没有因耗散作用而迅速消失, 而是在涡环内部维持着强度. 碰撞在充分发展后次级结构会从清晰的环状结构转变为模糊的湍流云. 构成次级结构的涡结构随着雷诺数的提升由分布清晰的涡丝变为紊乱的散点状涡区.
涡拟能同样受雷诺数影响显著. 研究发现, 随着雷诺数的提高, 初始时刻的涡拟能耗散速率呈现逐渐减小的趋势. 而由于碰撞过程涡核涡量提高, 涡拟能会呈现阶段性的爬升. 爬升的峰值随着雷诺数的提高而单调提高, 但是峰值对应的时间点则先延后再前移.
-
![]()
![]()
![]()
图 3 无扰动时涡环初始状态的特征量
Figure 3. Characteristic quantities of the initial state of the vortex ring in the absence of perturbations
![]()
图 4 涡结构变化结果对比$({{Re}} = 1600)$
Figure 4. Comparison between the vortex transformations $({{Re}} = 1600)$
![]()
图 5 动能发展结果对比(${Re} = 1000$)
Figure 5. Comparison between the kinetic energy developments (Re = 1000)
![]()
图 6 涡拟能发展结果对比$({Re} = 1000)$
Figure 6. Comparison between the enstrophy developments (Re = 1000)
![]()
图 7 不同扰动波数$ n $下涡量分布
Figure 7. Vorticity under different wave numbers $ n $ of perturbation
![]()
图 8 $T = 30$时雷诺数${Re} = 300$下涡量分布
Figure 8. Vorticity distribution under Reynolds number ${Re} = 300$ at $T = 30$
![]()
图 9 时间$T = 30$时碰撞中面涡量等值线
Figure 9. Collision midplane vorticity contours at time $T = 30$
![]()
图 11 时间$T = 8$时碰撞中面涡量等值线
Figure 11. Collision midplane vorticity contours at time $T = 8$
![]()
图 13 时间$T = 16$时涡量分布
Figure 13. Distribution of vorticity at different Reynolds numbers for $T = 16$
![]()
图 15 时间$T = 30$时碰撞中面涡量等值线
Figure 15. Collision midplane vorticity contours at $T = 30$
![]()
图 16 时间$T = 30$ ~ $ 37.5$涡量分布演变
Figure 16. Evolution of vortex distribution from time $T = 30$ ~ $37.5$
![]()
图 17 $T = 30$时不同雷诺数下涡量分布
Figure 17. Vorticity distribution at different Reynolds number for $T = 30$
![]()
图 18 T = 8不同雷诺数下碰撞中面涡量等值线
Figure 18. Collision midplane vorticity contours at different Reynolds number for T = 8
![]()
图 19 $T = 16$时不同雷诺数下涡量分布
Figure 19. Distribution of vorticity at different Reynolds number for $T = 16$
![]()
图 20 $T = 30$时不同雷诺数下涡量分布
Figure 20. Distribution of vorticity at different Reynolds number for $T = 30$
![]()
图 21 不同雷诺数下涡拟能在时间上的变化
Figure 21. Variation of normalized vortex fitting energy in time for different Reynolds numbers
![]()
图 22 涡拟能峰值高度随雷诺数的变化
Figure 22. The variation of peak enstrophy height with Reynolds number
![]()
图 23 涡拟能峰值对应时间随雷诺数的变化
Figure 23. The time position of peak enstrophy with Reynolds number
表 1 涡拟能峰值的独立性校验
Table 1 Independent verification of the peak value of enstrophy
dh dt 4.0 × 10−4 2.0 × 10−4 1.0 × 10−4 5.0 × 10−5 1/60 8.364 13.6964 8.6727 9.3375 1/80 13.6964 8.287 7.7339 8.9026 1/100 31.4801 8.5343 7.4551 8.2894 1/128 — 12.0731 6.5894 6.8445 1/150 — 19.575 5.3496 6.1625 
下载: 导出CSV
-
[1] 向阳, 刘洪, 吴镇远. 涡环物理特征的研究. 空气动力学学报, 2014, 32: 159-165 (Xiang Yang, Liu Hong, Wu Zhenyuan. The investigation on physical features of vortex rings. ACTA Aerodynamica Sinica, 2014, 32 (2): 159-165 (in Chinese) Xiang Yang, Liu Hong, Wu Zhenyuan. The investigation on physical features of vortex rings. ACTA Aerodynamica Sinica, 2014, 32 (2): 159-165 (in Chinese)
[2] 吴立新, 是勋刚. 理想流体中定常自由涡环运动. 力学学报, 1993, 25(5): 529-536 (Wu Lixin, Shi Xunguan. Steady free vortex rings in an inviscid fluid. Chinese Journal of Theoretical and Applied Mechanics, 1993, 25(5): 529-536 (in Chinese) Wu Lixin, Shi Xunguan. Steady free vortex rings in an inviscid fluid. Chinese Journal of Theoretical and Applied Mechanics, 1993, 25(5): 529-536 (in Chinese)
[3] 胡建军, 朱晴, 王美达等. 近距离下射流冲击平板PIV实验研究. 力学学报, 2020, 52(5): 1350-1361 (Hu Jianjun, Zhu Qing, Wang Meida, et al. PIV measurement of close impinging jet on flat plate. Chinese Journal of Theoretical and Applied Mechanics, 2020, 52(5): 1350-1361 (in Chinese) Hu Jianjun, Zhu Qing, Wang Meida, et al. PIV measurement of close impinging jet on flat plate. Chinese Journal of Theoretical and Applied Mechanics, 2020, 52(5): 1350-1361 (in Chinese)
[4] 徐杨, 王晋军. 涡环垂直撞击壁面的研究进展. 中国科学: 技术科学, 2013, 56: 2447-2455 (Xu Yang, Wang Jinjun. Recent development of vortex ring impinging onto the wall. Sci China Tech Sci, 2013, 56: 2447-2455 (in Chinese) doi: 10.1007/s11431-013-5338-7 Xu Yang, Wang Jinjun. Recent development of vortex ring impinging onto the wall. Sci China Tech Sci, 2013, 56: 2447-2455 (in Chinese) doi: 10.1007/s11431-013-5338-7
[5] Hernández RH, Reyes T. Symmetrical collision of multiple vortex rings. Physics of Fluids, 2017, 29(10): 103604 doi: 10.1063/1.5004587
[6] Tsai CY, Widnall SE. The stability of short waves on a straight vortex filament in a weak externally imposed strain field. Journal of Fluid Mechanics, 1976, 73(4): 721-733 doi: 10.1017/S0022112076001584
[7] Lim TT, Nickels TB. Instability and reconnection in the head-on collision of two vortex rings. Nature, 1992, 357(6375): 225-227 doi: 10.1038/357225a0
[8] Chu CC, Wang CT, Chang CC, et al. Head-on collision of two coaxial vortex rings: experiment and computation. Journal of Fluid Mechanics, 1995, 296: 39-71 doi: 10.1017/S0022112095002060
[9] Mckeown R, Ostilla-Mónico R, Pumir A, et al. Cascade leading to the emergence of small structures in vortex ring collisions. Physical Review Fluids, 2018, 3(12): 124702 doi: 10.1103/PhysRevFluids.3.124702
[10] Cheng M, Lou J, Lim TT. Numerical simulation of head-on collision of two coaxial vortex rings. Fluid Dynamics Research, 2018, 50(6): 065513 doi: 10.1088/1873-7005/aae54b
[11] Guan H, Wei ZJ, Rasolkova ER, et al. Numerical simulations of two coaxial vortex rings head-on collision. Advances in Applied Mathematics and Mechanics, 2016, 8(4): 616-647 doi: 10.4208/aamm.2014.m829
[12] Shariff K, Verzicco R, Orlandi P. A numerical study of three-dimensional vortex ring instabilities: viscous corrections and early nonlinear stage. Journal of Fluid Mechanics, 1994, 279: 351-375 doi: 10.1017/S0022112094003939
[13] Crow SC. Stability theory for a pair of trailing vortices. AIAA Journal, 1970, 8(12): 2172-2179 doi: 10.2514/3.6083
[14] Kerswell RR. Elliptical instability. Annual Review of Fluid Mechanics, 2002, 34(1): 83-113 doi: 10.1146/annurev.fluid.34.081701.171829
[15] Moore DW, Saffman PG. The instability of a straight vortex filament in a strain field. Proceedings of the Royal Society of London A Mathematical and Physical Sciences, 1975, 346: 413-425 doi: 10.1098/rspa.1975.0183
[16] Van Hooft JA, Popinet S. A fourth-order accurate adaptive solver for incompressible flow problems. Journal of Computational Physics, 2022, 462: 111251 doi: 10.1016/j.jcp.2022.111251
[17] Laporte F, Corjon A. Direct numerical simulations of the elliptic instability of a vortex pair. Physics of Fluids, 2000, 12(5): 1016-1031 doi: 10.1063/1.870357
[18] Cheng M, Lou J, Lim TT. Collision and reconnection of viscous elliptic vortex rings. Physics of Fluids, 2019, 31(6): 067107 doi: 10.1063/1.5095674
[19] Mishra A, Pumir A, Ostilla-Mónico R. Instability and disintegration of vortex rings during head-on collisions and wall interactions. Physical Review Fluids, 2021, 6(10): 104702 doi: 10.1103/PhysRevFluids.6.104702
[20] Melander M, Hussain F. Cut-and-connect of two antiparallel vortex tubes//Proceedings of the Summer Program, 1988: 257-286
[21] Yao J, Hussain F. Vortex reconnection and turbulence cascade. Annual Review of Fluid Mechanics, 2022, 54(1): 317-347 doi: 10.1146/annurev-fluid-030121-125143
[22] Cheng M, Lou J, Lim TT. Vortex ring with swirl: A numerical study. Physics of Fluids, 2010, 22(9): 097101 doi: 10.1063/1.3478976
[23] Hattori Y, Blanco-Rodríguez FJ, Le Dizès S. Numerical stability analysis of a vortex ring with swirl. Journal of Fluid Mechanics, 2019, 878: 5-36 doi: 10.1017/jfm.2019.621
[24] Moffatt HK. Generalised vortex rings with and without swirl. Fluid Dynamics Research, 1988, 3(1-4): 22 doi: 10.1016/0169-5983(88)90040-8
[25] Naitoh T, Okura N, Gotoh T, et al. On the evolution of vortex rings with swirl. Physics of Fluids, 2014, 26(6): 067101 doi: 10.1063/1.4882683
[26] Kollmann W, Prieto MI. DNS of the collision of co-axial vortex rings. Computers & Fluids, 2013, 73: 47-64
[27] Hussain F, Yao J. A physical model of turbulence cascade via vortex reconnection sequence and avalanche. Journal of Fluid Mechanics, 2020, 883: A51 doi: 10.1017/jfm.2019.905
[28] Rica S. Self-similar vortex reconnection. Comptes Rendus Mécanique, 2019, 347(4): 365-375
[29] Knio OM, Ghoniem AF. Numerical study of a three-dimensional vortex method. Journal of Computational Physics, 1990, 86(1): 75-106 doi: 10.1016/0021-9991(90)90092-F
[30] Mansfield JR, Knio OM, Meneveau C. Dynamic LES of colliding vortex rings using a 3D vortex method. Journal of Computational Physics, 1999, 152(1): 305-345 doi: 10.1006/jcph.1999.6258
计量
- 文章访问数: 258
- HTML全文浏览量: 44
- PDF下载量: 66
