Taylor-Couette模型是研究不稳定现象及流动转捩的经典模型[1]，因特殊的流动特性而被广泛应用，如黏度计、絮凝器、混合搅拌槽等．Taylor-Couette反应器由两个同心转筒构成，若保持外筒固定，则在内筒转速增加的过程中，两筒环隙间流体会依次经历层流涡、波状涡、调制波状涡和湍流涡等不同的流态转变[2]．文献[3]在研究中发现：波状涡流场可为絮凝反应提供最佳水力条件，但是在利用粒子图像速度场仪(PIV)[4]进一步挖掘研究波状涡流场的定量特征却遇到重重困难．不少学者对波状涡流态下的各种实验现象进行了研究．文献[5]对中性浮力粒子在Taylor-Couette流场三种流态下的惯性迁移进行了实验研究，结果表明波状涡流态下环隙内的粒子处于均匀混合的状态．文献[6]数值模拟了内外筒反向旋转的Taylor-Couette系统环隙间的流动情况，发现：在波状涡流态下，内外筒反向旋转与外筒固定内筒旋转的不同反应器内的波状涡流场的波速相似．文献[7]直接数值模拟了半径比为0.5的大间隙Taylor-Couette流动，发现旋转雷诺数(Re)为575时出现波状涡流动，且大间隙转捩过程中非轴对称扰动的出现很可能来源于外壁流动分离区域．这些属于宏观层面上的研究，对波状涡流场的变化细节还缺乏必要的研究．大涡模拟(the large eddy simulation，LES)结合直接数值模拟(DNS)与雷诺平均数值模拟(RANS)两者的优点，既保证了计算结果的精准性，也提高了计算效率[8]．不少学者采用大涡模拟方法[9-11]对Taylor-Couette流场开展了研究．本研究采用LES对Taylor-Couette反应器内的波状涡流场进行瞬态数值模拟，再现Taylor-Couette波状涡流场，以探究波状涡的变化规律及其周期变化特征．1 实验方法1.1　大涡模拟控制方程在求解大涡之前，要先过滤小于某一网格尺度的小涡，然后通过引入附加应力项来弥补小涡对大涡的影响[12]．对于不可压缩流动，经过滤波后的不可压缩纳维-斯托克斯方程即大涡模拟控制方程为      ∂u⃗i∂t+∂u⃗iu⃗j∂xj=-1ρ∂p⃗∂xi+ν∂2u⃗∂xjxj+∂∂xj(u⃗iu⃗j-uiuj⃗);∂u⃗i/∂xi=0；(u⃗iu⃗j-uiuj⃗)=τij，式中：ρ为流体密度；t为时间；u⃗i为i方向的速度矢量；u⃗j为j方向的速度矢量；uiuj⃗为流场速度分量乘积的矢量表达式；xi和xj为i，j方向的位移；u⃗为流场速度矢量；μ为流体黏性系数；p为流场的压强；τij为未知项，被称为亚格子雷诺应力，它是可解尺度脉动和不可解尺度脉动间的能量输运，须要采用Smagorinsky-Lilly亚格子模型进行封闭．亚格子应力定义为τij=τkkδij3-2μtSij，式中：τkk为亚格子湍动能；μt=(CsΔ)2|S|为涡黏系数，其中Cs为Smagorinsky常数，|S|=2SijSij为大尺度的应变值，△为网格尺度，Sij为经过滤波后的速度变形张量，又有Sij=12∂u⃗i∂xj+∂u⃗j∂xi．1.2　数值模型及数值方法Taylor-Couette模型几何尺寸如下：内筒半径r1=37.5 mm；外筒半径r0=50 mm；环隙宽度d=r0-r1=12.5 mm；内外筒半径比η=r1/r0=0.75；筒高L=440 mm．模型网格划分如图1所示，针对内外筒壁面，通过边界层理论计算得出y+值为1时的第一层网格高度y值，分别划分10层增长率为1.2的边界层网格．最后在满足网格无关性的要求下，综合考虑计算精度和计算速度，确定网格总数为2.09×106个．10.13245/j.hust.240816.F001图1Taylor-Couette模型网格划分示意图将划分好网格的Taylor-Couette模型文件导入到Fluent2020中进行相关设置．湍流模型选择LES，亚网格尺度模型选择Smagorinsky-Lilly，Cs取值在数值模拟方法验证后最终确定为0.1．边界条件如下：上下底面为自由表面；内外筒壁面分别为旋转和固定边界；压力速度耦合选择SIMPLE；空间离散化梯度项选择Least Square Cell Based；压力项为Second Order；动量项为Bounded Central Differencing；瞬态离散方案选择Second Order Implicit．1.3　涡识别方法与图像处理涡旋识别采用由文献[13]提出的Ω涡识别方法．该方法克服了传统涡识别方法须要针对具体流动调整阈值的问题，并能同时捕捉到强涡和弱涡．Ω涡识别方法将总涡量分解为旋转部分涡量R和非旋转部分涡量S，总涡量和旋转部分涡量的方向并不相同．引入参数Ω表示旋转部分涡量占总涡量的比例，Ω的表达式为Ω=BF2AF2+BF2+ε，式中：ε为一个小的正数，以避免当分母为极小数时出现误差极大的问题，文献[14]提出令ε=0.002Qmax为其近似值（Qmax为Q准则涡识别方法计算值最大值），很显然Ω的取值范围为0到1；A为对称张量；B为反对称张量．当Ω=1时，代表流体做刚体旋转；当Ω0.5时，代表B相较于A占优．因此可以采用Ω略大于0.5来作为涡识别的判据，Ω=0.8甚至Ω=0.9来寻找流场中的涡核结构．在数值模拟获得的环隙间流场的流动情况稳定后，截取环隙左侧子午面轴向坐标在190~250 mm范围内的涡作为研究对象，如图2所示，图中l为轴向距离．采用Image-Pro Plus 6.0对环隙子午面涡旋分布云图进行后处理，进而得到涡旋面积(图2中红色虚线内面积)及涡核的面积(图2中黑色虚线内面积)．10.13245/j.hust.240816.F002图2环隙子午面研究区域示意图1.4　工况选取Taylor-Couette流场流态的转变与旋转雷诺数Re相关[3]，有Re=ωr1d/ν，式中：ω为内筒旋转角速度；r1为内筒半径；d为环隙宽度；ν为流体的运动黏度(1.006×10-6 m2/s，20 ℃时水的运动黏度)．由于ω与内筒转速n的关系为ω=2πn，因此理论上控制内筒转速也能实现对旋转雷诺数的控制，旋转雷诺数的定义式可以改写为Re=πnr1d/(30ν)．由此可得：当Taylor-Couette反应器几何尺寸一定时，旋转雷诺数只与内筒转速有关．表1为工况及分析时间．根据前期的研究成果[3]，对转速在5~45 r/min (Re为244~2 196)范围内流场均进行了模拟．选取其中具有明显波状涡特征的流场进行分析，并给出部分工况及对应的分析时间．10.13245/j.hust.240816.T001表1工况及分析时间n/(r∙min-1)Re分析时间/s10，12，15，18488~878120~18020，25，28976~1 366120~15030，351 464~1 708120~140401 952120~1302 大涡模拟方法验证图3为不同旋转雷诺数下大涡模拟和PIV测量的环隙子午面速度矢量图对比，可以看出：大涡模拟得到的矢量图中涡旋形态与PIV测量结果相似，一定程度上可以反映大涡模拟结果的可靠性．而且大涡模拟结果比较规整有序，而PIV测量的矢量图可能是因为在测量过程中容易受到各种外界噪声因素干扰，显得略微杂乱，某些局部区域因噪声甚至会出现矢量紊乱的现象．综上比较看来：大涡模拟可以较为准确地反映出环隙间流场的流动特征，模拟结果具有一定的真实性和可靠性．10.13245/j.hust.240816.F003图3大涡模拟与PIV测量环隙子午面速度矢量图3 模拟结果根据模拟得到的不同时刻下的环隙子午面上的流场速度矢量场信息，采用Ω涡识别方法得到各时刻下环隙子午面上的涡旋、涡核分布云图．为在有限的篇幅内阐明问题，选取20，40 r/min的相关图示为代表，对转速为10 r/min(Re=486)，20 r/min(Re=972)，30 r/min(Re=1 459)，40 r/min(Re=1 945)的环隙子午面的涡旋分布云图及涡核分布云图进行分析，以研究波状涡的周期性变化规律．3.1　涡旋分布云图图4为Ω=0.52的环隙子午面上的涡旋分布云图周期展示．该图包含了各目标流场中一个完整的涡旋周期变化过程，通过5个时刻展示周期始、中、末的变化．选取各图中较为完整、清晰的一对涡旋作为目标，并分别编号为①号、②号涡旋进行观察．由系列图可以看到①号和②号涡旋经过一个完整的周期变化后又变回初始时刻的形状．10.13245/j.hust.240816.F004图4各工况环隙子午面涡旋分布云图周期变化3.1.1　涡旋大小观察①号、②号涡旋在系列云图中的变化过程，可以发现两相邻涡旋的面积大小发生了周期性的变化．当①号涡旋面积最大(②号涡旋面积最小)时，视作周期开始．从左至右，系列图中①号涡旋的面积由最大逐渐减小，同时②号涡旋的面积也由最小逐步增大；随着周期行进，①号涡旋的面积减至最小，②号涡旋的面积几乎同时增至最大；当周期结束时，①号涡旋的面积又增至最大，同时②号涡旋的面积也同步减至最小，两相邻涡旋经历一个周期变化后均变回周期开始时大小．当转速超出波状涡范围后，相邻涡旋大小几乎不再存在交替变化的现象．这种变化可以从能量变化的角度加以解释：当转速(Re)一定时，输入流场中的能量也是恒定的，涡旋大小的变化从一定程度上也反映了能量的动态分配与输运情况，所以根据能量守恒定律，流场中的涡旋大小要么保持恒定不变，要么以涡对中两相邻涡旋大小出现交替变化的互补形式以保持能量守恒．3.1.2　涡旋形态当转速为10和20 r/min时，随着周期行进，系列图中①号涡旋从最初的细长的芒果形逐渐收紧为粗短的芒果形，最后又拉伸回原先的细长的芒果形；与此同时，②号涡旋从最初的粗短的芒果形逐渐拉伸为细长的芒果形，最后又收紧回原先的粗短的芒果形．在涡旋拉伸的过程中，拉伸方向有趋近于指向旋转内筒的趋势，在云图中表现为向右上方或右下方偏移拉伸．当转速为30和40 r/min时，在一个周期内，两相邻涡旋的形态变化、拉伸幅度逐渐减小，变化差别逐渐减弱，最后的涡旋形状稳定地趋于橄榄形(甚至圆形)．高转速下各时刻的涡旋形态逐渐稳定并趋于相似，不再像较低转速时具有幅度很大的拉伸‐再收缩的形态变化过程．当转速超出波状涡范围后，相邻涡旋形态交替变化的现象难以捕捉到，几乎消失．3.1.3　涡心位置当转速为10和20 r/min时，随着周期行进，系列图中相邻涡旋的涡心位置在竖直方向和水平方向上发生了振动和漂移(竖直方向为振动，水平方向为漂移)，并在周期结束时刻又回到原点．当转速为30和40 r/min时，相邻涡旋的涡心位置在一个周期内虽也存在波动，但随着转速增加，其波动幅度越来越小，当转速超出波状涡范围后，涡心波动现象几乎消失．综上所述可知：当在三维空间中观察Taylor-Couette波状涡时，涡心位置会表现为前后左右上下的跳动；而在环隙子午面上，涡心位置随时间推移出现上下左右振动的现象．3.1.4　涡旋面积在分析取样时间内，每隔0.05 s采用Image-Pro Plus 6.0对各工况下的涡旋分布云图中的涡旋进行测量，①号和②号涡旋面积的测量结果见图5，图中：S1为涡旋面积；t为时间．10.13245/j.hust.240816.F005图5各工况涡旋面积周期变化图由图5可得：各目标流场中的涡旋面积存在周期性变化规律，对于每个目标流场，读取对应①号和②号涡旋面积变化曲线的峰值时间，将各个相邻峰值时间的极大值和极小值之间的差值取平均值，得到各目标流场中涡旋面积的变化周期．各工况的涡旋面积对应的变化周期分别为：10 r/min为12.94 s，20 r/min为6.80 s，30 r/min为1.93 s，40 r/min为1.49 s．可以看出：随着转速(n)增大，涡旋面积的变化周期减小．3.2　涡核分布云图图6为Ω=0.8的各工况下的环隙子午面上的涡核分布云图周期展示．该图包含了各目标流场中一个完整的涡核周期变化过程，通过5个时刻展示周期始、中、末的变化．选取各图中较为完整、清晰的一对涡核作为目标，并分别编号为①号、②号涡核．由系列图可以看到①号、②号涡核经过一个完整的周期变化后又变回初始时刻的形状．10.13245/j.hust.240816.F006图6各工况环隙子午面涡核分布云图周期变化过程3.2.1　涡核形态当转速为10 r/min时，涡核形态在一个周期内几乎不会发生变化，表现为斜向、小橄榄形．当转速为20 r/min时，随着周期行进，①号涡核经历由小变大再变小的过程，②号涡核则同步经历由大变小再变大的过程，整个周期变化过程中，涡核整体呈现斜向橄榄形，偶尔会出现涡核分散的不规则形状．当转速为30和40 r/min时，随着周期行进，①号涡核从最初的轴向、大橄榄形逐渐被挤压为径向、小橄榄形；与此同时，②号涡核从最初的径向、小橄榄形逐渐被拉伸为轴向、大橄榄形．可以看出：随着转速的增大，涡核面积也逐渐增大，涡核整体呈现出大而规则的形态，不再出现涡核分散等现象．3.2.2　涡核面积在分析取样时间内，每隔0.05 s采用Image-Pro Plus 6.0对各工况下涡核分布云图中的涡核进行测量，①号和②号涡核面积测量结果见图7，图中S2为涡核面积．10.13245/j.hust.240816.F007图7各工况涡核面积周期变化图由图7可知：各目标流场中的涡核面积存在周期性变化规律，对于每个目标流场，读取对应①号和②号涡核面积变化曲线的峰值时间，将各个相邻峰值时间的极大值和极小值之间的差值取平均值，得到各目标流场中涡核面积的变化周期．各工况的涡核面积对应的变化周期分别为：10 r/min为12.50 s，20 r/min为6.83 s，30 r/min为1.91 s，40 r/min为1.50 s．可以看出随着转速增大，涡核面积逐渐减小．4 分析与讨论4.1　波状涡周期与旋转雷诺数的关系为了探明涡旋与涡核的各相关周期特征指标之间及波状涡周期与内筒转速之间是否存在相关关系，将涡旋、涡核的周期与转速的关系作于图8，图中：T为各工况的周期时间；T1为涡旋周期时间；T2为涡核周期时间．10.13245/j.hust.240816.F008图8各工况下涡旋及涡核周期变化图图8为各工况下涡旋及涡核周期随内筒转速(Re)的变化情况．由图可知：各工况下波状涡变化周期随旋转雷诺数增大而逐渐减小；且无论是涡旋还是涡核，当Re一定时，二者周期几乎完全相同，这也说明涡核作为涡旋中强度最大的部分，直接影响了涡旋的整体性能特征，涡核的变化可以影响涡旋的变化，所以可以通过监测涡核的变化而预测涡旋的变化趋势．纵观Taylor涡的演变过程，在波状涡之后，随着旋转雷诺数增大，流场还要逐渐经历调制波状涡、湍流涡，最后完全变为湍流状态．在Re超出波状涡范围后，在调制波状涡区间内周期逐渐减小，直至最后难于捕捉；Re增大至湍流涡范围后，已经完全无法捕捉到周期现象．因此，周期随旋转雷诺数增大而减小，直至在湍流涡流场内周期消失的现象，可以说明周期时间是Taylor-Couette涡流场从有序向无序状态过渡的一种定量反映．有序而规律的变化随着旋转雷诺数的增大不断被压缩，体现在维持这种互补变化的周期逐渐减小，直至消失．最终，规律互补的变化也随着周期的消失而不复存在，涡旋完成由有序向无序状态的转变．所以周期的数值变化一定程度上可以反映并衡量涡流场由有序向无序状态的转变过程．4.2　涡核面积与涡旋面积之比λ涡核作为涡旋的一部分，为了探明二者之间是否存在某种联动或相互制约的关系，引入一个无量纲量λ．将λ定义为同一时刻涡核面积与涡旋面积之比，有λ=S1/S2．将一个周期内λ值随时间的变化情况作于图9，图中S为面积．另外，将涡旋、涡核的面积在一个周期中的变化情况同时在图中一并展示出来，以便于更直观地反映涡旋与涡核的同步变化关系，而且涡核在涡旋中的主导地位的变化过程也可以得到体现．10.13245/j.hust.240816.F009图9①号涡的λ与涡旋、涡核面积随时间的变化图9为各工况下①号涡的λ值与涡旋、涡核面积随时间的变化情况(②号涡变化趋势相反)．根据图9可知：在一个周期内，各工况下λ值与涡核面积随时间的变化趋势几乎一致，λ值主要受涡核面积的影响而与涡旋面积关系不大．说明λ值的变化趋势可代表涡核面积的变化趋势：随着转速的增加，λ值也增大，这意味着涡核面积同样也在变大，在涡旋中逐渐占据主导地位．观察图中涡旋与涡核的变化情况发现，在一个周期内，涡旋与涡核面积随时间的变化规律还与转速关系较大，不同转速下，两者的变化趋势不尽相同．当转速为10 r/min，20 r/min时，随着周期行进，涡旋与涡核面积的变化趋势并不一致，且几乎完全相反．这一规律从图4中也可以看出，即当涡核面积最大时，涡旋面积最小，且涡核面积随时间先增大后减小，涡旋面积却先减小后增大；在λ值上则表现为先增大后减小，这也同样反映出涡核与涡旋的面积大小呈现相反的变化规律．不能在表观几何层面上取得一致性说明存在能量方面因素的干扰：如果假设涡核是涡旋的能量中心，当涡核面积小时，能量中心比较弱，涡核凝聚力不足，涡旋展开发散，面积增大；当涡核面积大时，能量比较集中，能量中心强，涡核凝聚力强大，涡旋收拢集中，面积减小．涡旋与涡核前后时间段不一致的变化过程因为能量的引入而可以得到协调统一，即涡核带动涡旋在变化，两者是一个整体．整体性这一点在图8中也能得到体现，两者是合二为一的，体现在周期上几乎完全一致．所以，这种表观形貌在前后时间段变化不一致的现象可以认为是能量聚集与扩散互补的体现，λ的数值大小可以视作涡核能量强弱的一个反映．当转速为30 r/min，40 r/min时，随着周期行进，涡旋与涡核面积几乎同步增大或减小，变化趋势变得统一协调，只是相对来说，涡核增大或减小的幅度比涡旋更大一点．这与低转速下所表现出来的规律正好相反．这种现象同样可以引入能量因素加以解释：涡核是涡旋的能量中心，当旋转雷诺数增大后，涡核形态整体增大，几乎占据了整个涡区域，其能量更加集中，能量中心强度也不断增大，凝聚能力足够强大，对涡旋的整体把控也更加有力，不会再出现由于涡核面积减小、凝聚力不足而导致的涡旋表观面积扩散、涡旋大小虚增的现象．而在图中一个周期末位置处，涡核面积减小时刻涡旋面积仍在增大，可能是因为在过渡到下一个周期开始前，涡核已经先于涡旋发生变化，同样说明在周期变化中涡核处于主导地位．综合λ值与周期、涡旋、涡核之间的关系，均能说明涡核是涡旋的能量中心．在较低旋转雷诺数范围内，能量中心的强度相对较低，涡旋面积与涡核变化不同步，趋势相反：当涡核面积小时，能量中心弱，涡核凝聚力不足，涡旋面积大；当涡核面积大时，能量中心强，涡核凝聚力强，涡旋面积小．在较高旋转雷诺数下，能量中心的强度足够大，涡旋面积的变化完全受控于涡核的凝聚力，其大小变化同步依附于涡核．5 结论a．Taylor-Couette波状涡存在周期性变化规律，在一个周期中，环隙子午面上的涡旋表现出先增大后减小、或先减小后增大的变化过程．流场中存在涡对，相邻涡旋在一个周期中变化过程相反，一个膨胀变大，另一个收缩变小．b．Taylor-Couette波状涡流场中，涡旋、涡核面积的变化周期几乎完全相同，如当转速为20 r/min时，涡旋和涡核的周期变化时间分别为6.80和6.83 s；当转速为40 r/min时，涡旋和涡核的周期变化时间分别为1.49 s和1.50 s．涡旋、涡核两者共为一体，涡核作为涡旋中的强度核心在周期变化中处于主导地位，通过涡核的变化趋势可以预判涡旋的变化趋势．c．Taylor-Couette波状涡的周期随着旋转雷诺数的增大逐渐减小，说明周期的数值变化一定程度上可以反映并衡量Taylor-Couette涡流场形态的转变过程．