引 言
旋流分离器在石油、化工等行业有着广泛的应用．在井下油水分离过程中，与螺杆泵配合使用时，会使壁面绕其自身轴线旋转．旋流器主要分离段为变径圆管结构，内流场为强螺旋流，加之壁面的旋转，使流场更加复杂，是一种典型的强旋旋转湍流．目前，大量复杂的工程湍流问题仍采用基于Reynolds时均法的Reynolds应力模型进行计算[1]．相比涡黏模型，Reynolds应力模型由各向异性假设出发，能够描述更多的旋涡流动细节，同时又较大涡数值模拟或直接数值模拟等大大降低了资源及网格密度需求[2-3]．压力应变项的存在是Reynolds应力模型与涡黏模型最大的区别，也是不同Reynolds应力模型间的主要区别，是Reynolds应力模型的核心．但Sjogren等[4]和Mishra等[5]均指出，传统Reynolds应力模型中的压力应变项在应用于旋转湍流时，也会存在一定的问题，如不能给出正确的各向异性Reynolds应力张量的演化过程等，使得数值计算结果不够理想．Mishra等[5]提出只有通过在压力应变项模型中引入速度梯度的非线性相关项，才能解决Reynolds应力演化失真的问题．黄思源等[6]利用快速畸变理论对平面剪切与平面拉伸流动中Reynolds应力的演化特点，在压力应变快速项四阶张量的展开式中引入了平均尺度速度信息，对压力应变项进行了修正．
黄先北[7]在研究离心泵内旋转湍流时，也提出加入高阶非线性项能够使湍流模型更充分地体现旋转效应的影响．为了使Reynolds应力模型能够更好地描述壁面旋转变径圆管内的复杂旋转湍流问题，有必要对传统压力应变项进行适当修正．
1 壁面旋转变径圆管内流动特征分析
壁面旋转变径圆管内，螺旋流具有典型静态旋流器内流场的流动特点．同时，管壁绕其自身轴线的旋转又使其具备轴向旋转圆管内流体流动的部分特征．因此，由这两种基本的流动形式特征出发，对旋转条件下变径圆管内流场进行分析，为压力应变项的修正提供依据．
图1 变径圆管柱坐标系下速度的分解
Fig. 1 Velocity decomposition of the flow field in the reducing pipe in cylindrical coordinates
旋转变径圆管内螺旋流在柱坐标下的速度分解如图1所示，其中，vr，vz与vθ分别为径向速度、轴向速度与切向速度．螺旋流为强旋流动，切向速度最能反映出流场的流动特征，选取切向速度进行分析．
1.1 静态旋流器内流动的切向速度场特征
已有的研究表明，旋流器内切向速度场以最大切向速度为界分成内外两个涡流区，在外涡流区随半径减小切向速度逐渐增大，而在内涡流区随半径减小切向速度逐渐减小．
许多研究者拟合出水力旋流器内部的切向速度，确定了下面的关系[8]：
vθrn=const,
(1)
式中vθ为切向速度，r为半径，n为幂指数．当n=1时，为自由涡，在此区域内，旋转流体的流动为无旋流动，流体微团只沿曲线移动，而没有绕自身轴线的旋转运动，显然，此区域的流体切向速度vθ与半径r成反比．当n=-1时，为强制涡或似固旋转区，式(1)与固体旋转的表达式相同，强制涡可以认为是外涡流区流体连续高速旋转作用下造成的涡流，强制涡内流体微团不存在变形，只发生旋转．由上述分析可知，理论上旋流器的切向速度场是由外部的自由涡与内部的强制涡形成的组合涡．
根据切向速度场的实测结果，n并不是简单地等于1或者-1．由于实际流体的黏性作用，流体与壁面间的摩擦力以及切向速度沿径向的速度梯度过大的原因，使外涡流区的流体微团存在着一定的旋转，故外涡流区仅能称为准自由涡，但旋转作用在此区域并不占主导地位．内涡流区的流体微团也受到黏性等因素的影响，流体微团也会发生变形，其分布与强制涡有所不同，因此称为准强制涡．同样地，在准强制涡区，流体微团的旋转作用仍占主导地位．切向速度组合涡区结构[9]如图2所示．
1.2 轴向旋转圆管内流场特征
轴向旋转圆管典型的切向速度分布[10]见图3．即使当切向速度充分发展时，切向平均速度vθ分布，也不是强制涡结构，随着旋转数N的增大，vθ的速度剖面逐渐逼近抛物线，即
vθ=(r/R)2Vw,
(2)
式中Vw为管壁旋转速度，Vw=ωR．显然，圆管的旋转效应增强了切向方向上流体微团的旋转，与微团的变形效应相比，旋转效应明显占据了主导地位，流体的流动形态较未旋转的情况发生了改变．
图2 切向速度组合涡结构[9]图3 轴向旋转圆管切向速度[10]
Fig. 2 The composite vortex structure of Fig. 3 The tangential velocity in the the tangential velocity[9] axially rotating pipe[10]
当流动中流体微团的变形效应占主导地位，即变形率大于旋转率时，其流动称为双曲形流动；当微团的旋转效应占主导地位时，即旋转率大于变形率时，流动称为椭圆形流动；而当二者相等时，称为均匀剪切流动．从旋流器和绕轴旋转圆管的流场特征上看，旋转条件下变径圆管内部流场内涡区流动的旋转效应占主导地位，流动具有椭圆形流动的特征，外涡区的变形效应占主导地位，流动具有双曲形流动特征．
根据壁面旋转变径圆管内流体螺旋流动的分析，圆管轴截面上的速度分布具有椭圆流动的特征，在靠近轴心附近的内旋流区域内，速度分布甚至接近纯旋转流动，流体微团的旋转效应占据明显的主导地位．因此，将传统模型应用于此类型流动的模拟，将难以得到理想的结果．因此，需在现有压力应变项模型的基础上，发展一种能够适用于这类复杂旋转湍流模拟的压力应变项模型．
2 壁面旋转变径管内螺旋流的压力应变项修正模型
2.1 Reynolds应力模型
采用Reynolds时均法对壁面旋转变径管内螺旋流进行数值计算时，时均化后所产生的Reynolds应力需采用Reynolds应力模型进行计算．在均匀不可压湍流中，无量纲Reynolds应力的输运方程表示为
(3)
其中δij为Kronecker函数，Φij为压力应变项，目前已发展出多种形式的压力应变项模型，其中较为著名的有Launder等[11]推导出的线性(LRR)模型；Speziale等[12]提出的二阶应力(SSG)模型；Shih等[13]和Johansson等[14]提出的强可实现性的SL模型与JH模型等．
2.2 旋转湍流的压力应变项修正模型
针对旋转湍流问题，笔者已提出了相应的压力应变项修正模型，推导过程详见文献[15]．其基本思路是由快速畸变条件出发，分析传统压力应变项方程中缺少旋转率张量的相关量，无法描述旋转效应使流动趋于各向同性的问题．利用张量不变量理论，以衰减Reynolds应力各向异性张量为目的，构造出包含旋转率张量高阶的相关项．并以相同形式，构造出包含应变率张量高阶的相关项，将修正项由纯旋转流动拓展到微团旋转占优的旋转流动问题，建立了旋转湍流下Reynolds应力模型中压力应变项的修正模型，其形式如下：
(4)
式中
为旋转率张量高阶项；
为应变率张量高阶项；
为应变率张量；
为旋转率张量；
-BΩ=-ΩijΩji
为旋转率张量的不变量的相反数；
BS=SijSji
为应变率张量的不变量；
为各向异性Reynolds应力张量；
AbΩΩ=bikΩklΩlj, AbSS=bikSklSlj
均为混合不变量；
F=1-9Bb/2+9Cb
为湍流模型的可实现性条件约束；
Bb=bikbkl， Cb=bikbklbli
分别为各向异性Reynolds应力张量的第二与第三不变量．
2.3 针对壁面旋转变径圆管螺旋流的修正模型
由壁面旋转变径圆管内螺旋流的流动特征，其内涡流区为旋转占优的椭圆形流线流动，外涡流区为微团变形占优的双曲形流线流动．将上述修正模型应用与该流动时，应添加适当的模型参数，使其同时适用于内外涡流的模拟．
根据湍流模型的建模准则，模型系数应为不变量的函数．为了在系数中综合考虑流体微团旋转与变形对湍流Reynolds应力的影响，这里引入不变量β，其定义为
(5)
表示旋转率张量不变量在总不变量中所占的比例．对于微团应变主导的双曲流动，0<β<0.5，当流动为纯变形时，β=0，高阶项完全消失，不产生影响；而对于微团旋转主导的椭圆流动，0.5<β<1；而当纯旋转流动时，β=1．故根据流场中旋转与应变的大小动态变化，即有
(6)
添加包含速度梯度相关量高阶项的修正模型不仅可以正确描述旋转条件下湍流Reynolds应力的演化规律；同时由于管壁的旋转效应最直接影响到管内流体的速度梯度，模型中包含了速度梯度的高阶量，变径圆管壁面的旋转作用将通过影响周向速度的径向梯度形式进入旋转率张量与变形率张量，从而隐式地影响Reynolds应力的大小．尤其旋转率张量包含了壁面所诱导产生的沿轴向方向上的涡量分量，其会直接影响Reynolds应力．
3 修正模型用于旋转变径圆管螺旋流流场的验证计算
3.1 数学模型
以SSG模型为基础，其表达式如下：
(7)
式中
为模型系数．将SSG模型与修正模型联立，得到修正后的压力应变项如下：
(8)
3.2 几何模型
以56 mm Thew型旋流器为变径圆管基本结构，几何模型结构如图4所示(D=56 mm，Dl/D=0.5,Du/D=0.25,Do/D=1/14,Lc/D=1.6,Lu/D=5,α=20°,β=1.5°)．流体从双侧切向入口流入，形成螺旋流，经过旋流腔、大锥段、小锥段与尾管段，内涡流区流体从溢流出口流出，外涡流区流体从底流出口流出．旋转条件下，变径圆管按图示方向整体绕自身轴线旋转．
图4 几何模型
Fig. 4 The geometric model
3.3 计算网格与边界
Reynolds应力模型对近壁区的求解需借助壁面函数法，按壁面函数法的求解要求，把第一层网格中心布置在对数率层，使y+≈15，并保证边界层内网格有足够的分辨率．网格总数为199×64×23(轴向×周向×径向)．计算域除中心少量楔形体网格外，均为六面体结构网格．
以N定义壁面旋转条件，N=Vw/Um为无量纲旋转率．以主直径D处Reynolds数定义流动条件．定义溢流分流比，具体参数如表1所示．
表1 计算边界条件
Table 1 Boundary conditions
descriptionsymbolvaluemean velocity at DUm/(m/s)0.325tangential velocity of wall at DVw/(m/s)1.46Reynolds number at DRe19 000rotation rateN4.5split ratioε/%5
3.4 离散格式与算法
控制方程在空间上采用有限体积法进行离散，具体离散格式及算法见表2．
表2 离散格式与算法
Table 2 The discrete scheme and algorithm
discrete scheme or algorithm methodgradient algorithmleast squares methodpressure interpolation algorithmstandardconvection interpolation2nd-order upwindpressure-velocity coupling methodSIMPLE
3.5 计算结果与分析
图5与图6分别给出了z=150，200，300 mm共3个截面上修正模型与未修正模型计算的切向速度、轴向速度与实测结果的对比，截面位置如图4所示．实测结果来源于激光Doppler测速(LDV)方法，具体测试方法详见文献[16]．
(a) z=150 mm(b) z=200 mm(c) z=300 mm
图5 切向速度对比
Fig. 5 Comparison of tangential velocities
可以看出，3种方法得到切向速度分布趋势基本上是一致的，但在数值上，两种模型的切向速度在数值上都偏小．相比SSG，COR的最大切向速度点更靠外．在外涡流区，COR计算结果的数值更大，数值上与实测值更接近；内涡流区，SSG计算结果数值更大，更接近实测．分析计算结果与实测存在差异的原因，主要有以下两点：首先，修正模型仍基于Reynolds时均法与Reynolds应力模型，计算精度受到了方法本身精度的限制；其次，已有文献表明，壁面的旋转效应会改变近壁区流动的边界层结构[17]，而现有传统壁面函数对其模拟能力有限，本文采用传统的壁面函数对近壁区流动进行描述，因此，造成外涡流区的切向速度与实测仍有一定的差距．但总体上看，修正模型的切向速度在整体趋势上与实测结果更接近．
(a) z=150 mm(b) z=200 mm(c) z=300 mm
图6 轴向速度对比
Fig. 6 Comparison of axial velocities
不同模型轴向速度的数值计算结果与实测结果在数值上都较为接近．两种模型在内外涡区的计算结果同样出现了一定的差异，在外涡流区，SSG与COR在数值上差别不大，其分布规律与实测值基本上也是一致的；而在内涡流区，SSG与COR产生了明显的差异，COR的计算结果无论从数值上还是趋势上，都明显与实测值更为接近．
4 结 论
1) 根据旋流器与轴向旋转圆管内流场特征，对旋转变径圆管内螺旋流的特征进行了分析，确定了压力应变项修正的依据．
2) 基于适用于微团旋转占主导流动的修正压力应变项，引入了不变量形式的模型系数，综合考虑流体微团旋转与变形对湍流Reynolds应力的影响，且随着微团旋转与变形间的大小关系而变化．
3) 将修正项模型与二阶应力压力应变项模型联立，构建了新的压力应变项．采用修正后模型模拟研究了壁面旋转变径圆管内螺旋流动，将未修正、修正后及实测切向速度、轴向速度进行了对比，修正模型在速度场的预测上要优于未修正模型，其结果更接近实测结果．对于这类复杂旋转湍流而言，非线性的修正模型对传统模型起到了一定的改进作用．
[1] 冯志鹏, 臧峰刚, 张毅雄. 三维横向流体诱发直管振动的数值模拟[J]. 应用数学和力学, 2013, 34(12): 1311-1320.(FENG Zhipeng, ZANG Fenggang, ZHANG Yixiong. Numerical simulation of uniform tube vibration induced by 3-dimensional cross-flow[J]. Applied Mathematics and Mechanics, 2013, 34(12): 1311-1320.(in Chinese))
[2] 王福军. 计算流体动力学分析[M]. 北京: 清华大学出版社, 2004.(WANG Fujun. Computational Fluid Dynamics Analysis[M]. Beijing: Tsinghua University Press, 2004.(in Chinese))
[3] XU C Y, ZHOU T, WANG C L, et al. Applications of scale-adaptive simulation technique based on one-equation turbulence model[J]. Applied Mathematics and Mechanics (English Edition), 2015, 36(1): 121-130.
[4] SJOGREN T, JOHANSSON A V. Development and calibration of nonlinear algebraic models for terms in Reynolds stress transport equations[J]. Physics of Fluids, 2000, 12(6): 1554-1572.
[5] MISHRA A A, GIRIMAJI S S. Pressure-strain correlation modeling: towards achieving consistency with rapid distortion theory[J]. Flow Turbulence Combustion, 2010, 85(3/4): 593-619.
[6] 黄思源, 符松. 快速畸变下压力应变快速项模式研究[J]. 中国科学 (G辑), 2008, 38(9): 1255-1264.(HUANG Siyuan, FU Song. Study on the rapid strain and pressure term model under rapid distortion[J]. Science in China (Series G), 2008, 38(9): 1255-1264.(in Chinese))
[7] 黄先北. 动态三阶非线性SGS模型及其在离心泵流动模拟中的应用[D]. 博士学位论文. 北京: 中国农业大学, 2017.(HUANG Xianbei. Dynamic cubic non-linear SGS model and the application in centrifugal pumps[D]. PhD Thesis. Beijing: China Agricultural University, 2017.(in Chinese))
[8] BRADLEY D, PULLING D J. Flow patterns in the hydraulic cyclone and their interpretation in terms of performance[J]. Transactions Institution of Chemical Engineers, 1959, 37: 34-45.
[9] 蒋明虎, 赵立新, 李枫, 等. 旋流分离技术[M]. 哈尔滨: 哈尔滨工业大学出版社, 2000.(JIANG Minghu, ZHAO Lixin, LI Feng, et al. Hydrocyclone Separation Technology[M]. Harbin: Harbin Institute of Technology Press, 2000.(in Chinese))
[10] IMAO S, ITOH M, HARADA T. Turbulent characteristics of the flow in an axially rotating pipe[J]. International Journal of Heat and Fluid Flow, 1996, 17(5): 444-451.
[11] LAUNDER B E, REECE G J, RODI W. Progress in the development of a Reynolds-stress turbulence closure[J]. Journal of Fluid Mechanics, 1975, 68(3): 537-566.
[12] SPEZIALE C G, SARKAR S, GATSKI T B. Modelling the pressure-strain correlation of turbulence: an invariant dynamical systems approach[J]. Journal of Fluid Mechanics, 1991, 227: 245-272.
[13] SHIH T H, LUMLEY J L. Modeling of Pressure Correlation Terms in Reynolds Stress and Scalar Flux Equations[M]. Cornell University, 1985.
[14] JOHANSSON A V, HALLB
CK M, LINDBORG E. Modelling of rapid pressure-strain in Reynolds stress closures: difficulties associated with rotational mean flows[J]. Applied Scientific Research, 1994, 53(1/2): 119-137.
[15] 张井龙, 王尊策, 徐艳. 旋转湍流中雷诺应力模型压力应变关联项的修正研究[J]. 计算力学学报, 2014, 31(6): 787-792.(ZHANG Jinglong, WANG Zunce, XU Yan. Pressure strain correlation term of Reynolds stress model for rotating turbulence[J]. Chinese Journal of Computational Mechanics, 2014, 31(6): 787-792.(in Chinese))
[16] 徐艳, 张彦月, 张井龙, 等. 基于LES方法的壁面旋转旋流分离器内流特性的研究[J]. 流体机械, 2017, 45(8): 1-5.(XU Yan, ZHANG Yanyue, ZHANG Jinglong, et al. Study on the internal flow characteristics of a wall rotating hydrocyclone based on LES method[J]. Fluid Machinery, 2017, 45(8): 1-5.(in Chinese))
[17] 刘难生. 旋转湍流的直接数值模拟和大涡模拟研究[D]. 博士学位论文. 合肥: 中国科学技术大学, 2003: 88-89.(LIU Nansheng. Direct numerical simulation and large eddy numerical simulation of rotating turbulence[D]. PhD Thesis. Hefei: University of Science and Technology of China, 2003: 88-89.(in Chinese))