崔振东^1,2, 刘斌¹, 蔡军¹, 淮秀兰¹
1. 中国科学院工程热物理研究所, 北京 100190;
2. 中国科学院大学, 北京 100190
2015年03月13日 收稿; 2015年05月25日 收修改稿
基金项目: 国家自然科学基金(51376181)资助
通信作者: E-mail: caijun@mail.etp.ac.cn

摘要: 应用4种湍流模型(标准k-ε模型、标准k-ω模型、Transition SST模型和Transition k-kl-ω模型)对微通道内空化流动进行模拟.通过与实验结果对比,分析4种湍流模型在空化流动模拟计算中的预测能力.结果表明:标准k-ω模型预测的流量和空化流型与实验结果最为接近.Transition SST模型和Transition k-kl-ω模型流量预测结果与实验较为吻合,但流型预测误差较大.标准k-ε模型未能预测出空化流动,不适用于低雷诺数下的空化流动计算.
关键词: 水力空化微通道湍流模型数值模拟
Analysis of simulation results for cavitating flows in micro-channels using different turbulence models
CUI Zhendong^1,2, LIU Bin¹, CAI Jun¹, HUAI Xiulan¹
1. Institute of Engineering Thermophysics, Chinese Academy of Sciences, Beijing 100190, China;
2. University of Chinese Academy of Sciences, Beijing 100190, China


Abstract: The cavitation flows in micro-channels were numerically investigated using four different turbulence models, and the models are as follows: standard k-ε turbulence model, standard k-ω turbulence model, Transition SST turbulence model, and Transition k-kl-ω turbulence model. In comparison to the experimental results, the prediction abilities of these turbulence models for the cavitation flow were analyzed in details. The results show that both the flow rate and flow pattern predicted by the standard k-ε turbulence model agree well with the experimental results. Although the flow rates predicted by the Transition SST turbulence model and the Transition k-kl-ω turbulence model also agree well with the experimental results, the predicted flow patterns significantly deviate from the experimental ones. The standard k-ε turbulence model failed to predict the cavitation flow, which indicates that it is not suitable for the low Reynolds number flow circumstance.
Key words: hydrodynamic cavitationmicro-channelturbulence modelnumerical simulation
当液体流经一定的水力限流结构，如孔板和文丘里管等，会在过流表面局部产生低压，从而诱发空化^[1].在流体设备中，空化泡的溃灭会产生高压脉冲及微射流，从而造成壁面材料的侵蚀^[2].近年来，随着微加工技术的日趋成熟，微流体系统广泛应用于航空航天、生物医药和仪器检测等领域^[3-4].与此同时，如何消除空化作用对微流体系统工作的影响，成为研究热点.各国****通过实验方法对微尺度下的空化流动特性进行了广泛研究^[5-9].但是，由于空化现象是一种复杂的两相流动，当前的实验手段很难获得详细的流动参数.因此，对微尺度下空化流动的认识更多地需要依靠数值模拟.Hickel等^[10]首次采用大涡模拟(LES)的方法对微通道内的空化流动进行模拟，所得结果与实验吻合良好.但由于运算量巨大，LES只能应用于简单空化流动的预测.同样，受计算速度和容量限制，直接数值模拟(DNS)方法也无法实际应用.目前的数值研究主要基于Reynolds平均数值模拟(RANS)方法.Rooze等^[11]应用SST湍流模型对内嵌矩形限流结构微通道空化流动进行模拟，计算所得流型与实验取得了定性上的一致.然而，湍流模型用于微通道内空化流动的研究还处于起步阶段，计算精度有待提高.本文采用常用的4个湍流模型，对微通道内空化流动进行数值研究，比较分析不同湍流模型对流动预测的适应性和准确性，以期为微通道空化流动模拟湍流模型的选择和计算结果的理解提供理论依据.
1 物理及数学模型1.1 物理模型选取内嵌限流结构的矩形微通道为研究对象，结构形式如图 1(a)所示.根据微通道在z方向上的对称性，选取通道结构的一半进行计算，以提高计算效率.以去离子水作为工作介质，在进、出口压力差的作用下，依次流经微通道入口段、限流结构段和微通道出口段.物理模型计算参数：微通道高度H=0.2mm，微通道半宽度W=0.25mm，微通道入口段长度L₁=2mm，限流结构长度L₂=0.2mm，微通道出口段长度L₃=23.8mm.模型采用结构化网格划分方法，同时对限流结构局部网格进行加密处理，以保证对流场的精确捕捉，如图 1(b)所示.为消除网格数目对计算结果的影响，对网格无关性进行了分析.结果显示，当网格数大于1.48×10⁵时，模拟结果变化不大，可以认为满足要求.
Fig. 1

	Download: JPG larger image
图 1 内嵌限流结构的矩形微通道结构示意图 Fig. 1 Schematic of the micro-channel with rectangular orifice

1.2 数学模型稳态条件下，空化流动过程中的质量、动量和能量守恒方程分别表述为
1) 质量守恒方程

${{\partial {\rho _m}} \over {\partial t}} + {\partial \over {\partial {x_j}}}({\rho _m}{u_j}) = 0,$

(1)

式中，t为时间，u_j为速度，ρ_m为气液混合物密度，根据流体体积百分比进行计算

${\rho _m} = {\rho _1}{\alpha _1} + {\rho _v}{\alpha _v},$

(2)

其中，α为体积分数，下标l和v分别代表液相和气相.
2) 动量守恒方程

${\partial \over {\partial t}}({\rho _m}{u_i}) + {\partial \over {\partial {x_j}}}({\rho _m}{u_i}{u_j}) = - {{\partial p} \over {\partial {x_i}}} + {\rho _m}{g_i} + {{\partial {\tau _{ij}}} \over {\partial {x_i}}}$

(3)

式中，μ_m为混合物动力粘性系数，τ_ij为剪切应力，表达式为

${\tau _{ij}} = \mu \left( {{{\partial {u_i}} \over {\partial {x_j}}} + {{\partial {u_j}} \over {\partial {x_i}}} - {2 \over 3}{\delta _{ij}}{{\partial {u_k}} \over {\partial {x_k}}}} \right).$

(4)

3) 空化模型
本文采用Schnerr和Sauer发展的空化模型^[12]，蒸汽输运方程表示为

${\partial \over {\partial t}}({\rho _v}{\alpha _v}) + {\partial \over {\partial {x_i}}}({\rho _v}{u_i}{\alpha _v}) = R$

(5)

式中，α_v为蒸汽相体积分数，R为蒸汽净生成速率，其表达式为

$R = {{{\rho _v}{\rho _1}} \over {{\rho _m}}}\left[ {{{\partial {\alpha _v}} \over {\partial t}} + {{\partial ({\mu _j}{\alpha _v})} \over {\partial {x_j}}}} \right].$

(6)

空泡半径R_b为

$Rb = {{{\alpha _v}} \over {1 - {\alpha _v}}}{3 \over {4/\pi }}{1 \over {{n_b}}},$

(7)

式中，n_b是模型常数，本文取值为1.2×10¹⁰，因此可得

$R = {{{\rho _v}{\rho _1}} \over {{\rho _m}}}{\alpha _v}(1 - {\alpha _v}){3 \over {{R_b}}}{\rm{sign}}({p_{sat}} - p)\sqrt {{2 \over 3}{{\left| {{p_{sat}} - p} \right|} \over {{\rho _1}}}} ,$

(8)

其中，p_sat为饱和蒸汽压.
4) 湍流模型
标准k-ε湍流模型中，湍动能k和湍动耗散率ε的输运方程如下

$\eqalign{ & {\partial \over {\partial t}}({\rho _m}k) + {\partial \over {\partial {x_i}}}({\rho _m}k{u_i}) = \left[ {{\partial \over {\partial {x_j}}}({\mu _m} + {{{\mu _t}} \over {{\sigma _k}}}){{\partial k} \over {\partial {x_j}}}} \right] + \cr & {G_k} + {G_b} - {\rho _m}\varepsilon - {Y_M} + {S_k}, \cr} $

(9)

$\eqalign{ & {\partial \over {\partial t}}({\rho _m}\varepsilon ) + {\partial \over {\partial {x_i}}}({\rho _m}\varepsilon {u_i}) = {\partial \over {\partial {x_j}}}\left[ {({\mu _m} + {{{\mu _t}} \over {{\sigma _\varepsilon }}}){{\partial \varepsilon } \over {\partial {x_j}}}} \right] + \cr & {C_1}_\varepsilon \varepsilon k({G_k} + {C_1}_\varepsilon Gb) - {C_{2\varepsilon }}\rho {{{\varepsilon ^2}} \over k} + {S_\varepsilon }, \cr} $

(10)

式中，G_k为速度梯度引起的湍动能，表达式为${G_k} = {\mu _i}\left( {{{\partial {\mu _i}} \over {\partial {x_j}}} + {{\partial {\mu _j}} \over {\partial {x_i}}}} \right){{\partial {\mu _i}} \over {\partial {x_j}}};Gb$为浮力引起的湍动能，表达式为$Gb = \beta {g_i}{{{\mu _t}} \over {P{r_t}}}{{\partial T} \over {\partial {x_i}}};\beta $是热膨胀系数；g_i是重力加速度在第i方向的分量；Pr_t是湍动普朗特数；Y_M为可压缩湍流中脉动扩散产生的波动，表达式为Y_M=2ρεM_t²；M_t是湍流马赫数；S_k和S_ε为用户定义的源项，本文中取值为0；μ_t为湍流黏度.
标准k-ω湍流模型中，湍动能k方程和耗散率ω的输运方程如下

$\eqalign{ & {\partial \over {\partial t}}({\rho _m}k) + {\partial \over {\partial {x_i}}}({\rho _m}k{u_i}) = {\partial \over {\partial {x_j}}}\left[ {{\Gamma _k}{{\partial k} \over {\partial {x_j}}}} \right] + \cr & {G_k} - {Y_k} + {S_k}, \cr} $

(11)

$\eqalign{ & {\partial \over {\partial t}}({\rho _m}\omega ) + {\partial \over {\partial {x_i}}}({\rho _m}\omega {u_i}) = {\partial \over {\partial {x_j}}}\left[ {{\Gamma _\omega }{{\partial \varepsilon } \over {\partial {x_j}}}} \right] + \cr & {G_\omega } - {Y_\omega } + {S_\omega }, \cr} $

(12)

式中，G_ω为湍流耗散率的生成项，Г_k和Г_ω分别为k和ω的有效扩散率，Y_k和Y_ω分别为k和ω的耗散.上述变量的具体表达式参见文献^[13].
Transition SST湍流模型是在SST k-ω模型的基础上增加了对湍动能k的修正，并且与其他2个输运方程耦合而产生的.
间歇因子γ的输运方程为

$\eqalign{ & {\partial \over {\partial t}}({\rho _m}\gamma ) + {\partial \over {\partial {x_j}}}({\rho _m}{U_j}\gamma ) = {P_{\gamma 1}} - {E_{\gamma 1}} + {P_{\gamma 2}} - \cr & {E_{\gamma 2}} + {\partial \over {\partial {x_j}}}\left[ {\left( {\mu + {{{\mu _t}} \over {{\sigma _\gamma }}}} \right){{\partial \gamma } \over {\partial {x_j}}}} \right]. \cr} $

(13)

转捩开始的动量厚度雷诺数Re_θt的输运方程为

$\eqalign{ & {\partial \over {\partial t}}({\rho _m}{\mathop {Re}\limits^ \sim _{\theta t}}) + {\partial \over {\partial {x_j}}}({\rho _m}{U_j}{\mathop {Re}\limits^ \sim _{\theta t}}) = {P_{\theta t}} + \cr & {\partial \over {\partial {x_j}}}\left[ {{\sigma _{\theta t}}(\mu + {\mu _t}){{\partial {{\mathop {Re}\limits^ \sim }_{\theta t}}} \over {\partial {x_j}}}} \right]. \cr} $

(14)

修正后的湍动能k方程如下

$\eqalign{ & {\partial \over {\partial t}}({\rho _m}k) + {\partial \over {\partial {x_i}}}({\rho _m}k{u_i}) = {\partial \over {\partial {x_j}}}\left( {{\Gamma _k}{{\partial k} \over {\partial {x_j}}}} \right) + \cr & {G^*}_k + {Y^*}_k + {S_k}, \cr} $

(15)

式中，G_k^*和Y_k^*的计算方法为

${G^*}_k = {\gamma _{{\rm{eff}}}}{{\tilde G}_k}$

(16)

${Y^*}_k = {\rm{min}}({\rm{max}}({\gamma _{{\rm{eff}}}},0.1){\rm{ }},1.0){\rm{ }}{Y_k}.$

(17)

Transition k-kl-ω转捩模型是用来预测边界层发展的湍流模型，包含湍流动能k_T、层流动能k_L和耗散率ω输运方程

$\eqalign{ & {{D{k_T}} \over {Dt}} = {P_{KT}} + R + {R_{NAT}} - \omega {k_T} - {D_T} + \cr & {\partial \over {\partial {x_j}}}\left[ {\left( {v + {{{\alpha _T}} \over {{\alpha _k}}}} \right){{\partial {K_T}} \over {\partial {x_j}}}} \right], \cr} $

(18)

${{D{k_T}} \over {Dt}} = {P_{KL}} - R - {R_{NAT}} - {D_L} + {\partial \over {\partial {x_j}}}\left[ {v{{\partial {K_L}} \over {\partial {x_j}}}} \right]$

(19)

$\eqalign{ & {{D\omega } \over {Dt}} = {C_{\omega 1}}{\omega \over {{k_T}}}{P_{kT}} + \left( {{C_{\omega R}}{f_\omega } - 1} \right){\omega \over {{k_T}}}\left( {R + {R_{NAT}}} \right) - \cr & {C_{\omega 2}}{\omega ^2} + {C_{\omega 3}}{f_\omega }{\alpha _T}f_W^2{{\sqrt {kT} } \over {{d^3}}} + {\partial \over {\partial {x_j}}}\left[ {\left( {v + {{{\alpha _T}} \over {{\alpha _k}}}} \right){{\partial \omega } \over {\partial {x_j}}}} \right], \cr} $

(20)

式中，P_KT为湍流动能产生源项，P_KL为大尺度湍流脉动引起的层流动能源项，R为由层流动能所转化的湍动能的生成率，R_NAT为自然转捩的产生源项，上述参数的计算表达式与孙润鹏等^[14]的相同.上述各湍流模型中计算常数均选取默认值.
2 数值方法空间离散方法采用有限体积法，应用耦合求解器对控制方程进行求解.压力离散采用PRESTO格式，动量方程、蒸汽输运方程与湍流方程均采用绝对稳定的一阶迎风格式.液相为去离子水，气相为饱和蒸汽.进、出口边界均采用压力边界条件.固壁边界采用无滑移速度边界，黏性底层采用标准壁面函数进行处理.
3 结果与讨论为便于分析，首先定义表征空化流动特性的无量纲参数，即雷诺数(Re)和空化数(σ)，其表达式为

$Re = {{\rho \bar uDh} \over \mu }$

(21)

$\sigma = {{{p_{{\rm{out}}}} - {p_v}} \over {{1 \over 2}\rho {{\bar u}^2}}}$

(22)

式中，D_h为微通道限流结构当量直径.
当空化数低于临界值时，通道内流体将发生空化现象，并且随着空化数的降低，空化现象更为明显.本文参照前期实验数据^[15]，选取雷诺数Re=2405，空化数为σ=0.0348的工况进行模拟，对应进、出口压力参数如表 1所示.
Table 1

表 1 计算参数Table 1 Calculation parameters 空化数 pin /kPa pout /kPa σ=0.0348 355 8.5

表 1 计算参数Table 1 Calculation parameters

4种湍流模型计算所得的流量及其与实验数据的对比如表 2所示.可以看出，标准k-ω湍流模型流量预测结果与实验数据最为接近，相对误差约为4.59%；Transition SST和Transition k-kl-ω湍流模型所得流量与实验结果也较为吻合，相对误差分别为5.09%和5.05%；标准k-ε湍流模型所得流量与实验值偏差最大，相对误差达到8.75%.由于标准k-ε湍流模型基于完全湍流假设，只适合于湍流充分发展的流动过程模拟，因此对本算例低雷诺数空化流动的预测精度较差.标准的k-ω模型是一种低雷诺数模型，考虑了剪切流效应，如尾流、绕流、和放射状喷射流等，对墙壁束缚流动和自由剪切流动模拟效果较好.微通道空化流动中，在限流结构下游出现的尾焰状气-液两相射流受到壁面的显著影响，因此采用标准k-ω模型计算的流量具有较高的精度.
图 2所示为微通道y=0.1mm截面处4种湍流模型所得空化流型与实验结果的对比.图 2(a)为实验观察到的空化充分发展时的流动形态.可以看出，在限流结构下游，微通道中轴线附近为液体射流，其两侧是稳定的空化泡，里面充满饱和水蒸气.在空化区域下游，由于压力的恢复，空化泡的溃灭导致汽含率迅速降低，下游通道内的流动最终恢复单相.标准k-ε湍流模型计算所得汽含率始终为0，未能预测出空化现象，见图 2(b)所示.结合表 2流量结果可知，对于目前微通道结构，k-ε湍流模型不适用于Re≤2405的空化流动计算.图 2(c)所示为k-ω湍流模型所得汽含率分布.由于空化区域气-液两相的存在，流体中的分子粘性对流态产生重要作用，而k-ω湍流模型能准确的捕捉到壁面附近黏性底层的剪切流动，因此所得空化流型与实验结果吻合良好.Transition SST和Transition k-kl-ω湍流模型所得的汽含率分布分别如图 2(d)和2(e)所示.2个模型的共同点是考虑了转捩机理，但是模拟所得的空化射流过长，与实验结果存在较大差距.由此可知，对于充分发展的空化流动，整个流场中湍流占据主导地位，无需考虑层流向湍流转变的影响.
Table 2

表 2 不同湍流模型计算所得流量与实验结果对比Table 2 Comparison of flow rate between the experimental results and the calculated results using four turbulence models 项目 流量 相对误差/% 实验数据 9.60 —— k-ε 8.76 8.75 k-ω 9.16 4.59 Transition SST 9.11 5.09 Transition k-kl-ω 9.12 5.05

表 2 不同湍流模型计算所得流量与实验结果对比Table 2 Comparison of flow rate between the experimental results and the calculated results using four turbulence models

Fig. 2

	Download: JPG larger image
图 2 微通道y=0.1 mm截面空化流场汽含率分布 Fig. 2 Gas fraction distributions of the cavitation flow field in the micro-channel taken in the plane at y=0.1mm

图 3所示为微通道y=0.1mm截面处不同湍流模型所得速度场的分布.对比不同模型结果可以发现，在限流结构下游两侧壁面附近均出现回流区，其中k-ω湍流模型预测的回流区面积和流速均明显高于其他3个模型.结合图 2可以看出，回流区内的汽含率最高，与空化发生位置相吻合，二者之间具有紧密联系.由此可以判断，回流区内流动预测的正确与否对空化流场计算结果的准确性有至关重要的影响.由于k-ω湍流模型在回流计算方面相比其他3个模型具有更好的适应性，因此模拟结果与实验最为接近.
Fig. 3

	Download: JPG larger image
图 3 微通道y=0.1 mm截面速度场分布 Fig. 3 Velocity distributions in the micro-channel taken in the plane at y=0.1mm

4 结论本文将标准k-ε模型、标准k-ω模型、Transition SST模型和Transition k-kl-ω模型用于微通道内空化流动模拟，并将计算结果与实验数据进行对比，得出如下结论：
1) 标准k-ω模型、Transition SST和Transition k-kl-ω湍流模型预测的流量结果均较为接近实验数据，其中标准k-ω模型吻合最好，而标准k-ε模型所得结果与实验数据存在较大偏差.
2) 气-液相变中流体粘性对流型产生重要影响，标准k-ω模型能较好预测壁面附近黏性底层的剪切流动，所得流型结果与实验最为接近，其他3个模型结果均存在较大差距.
3) 微通道限流结构下游两侧壁面附近存在回流区，与空化发生位置相吻合，模型对回流区预测的准确性影响到计算结果的精度.
参考文献

[1]	王献孚. 空化泡和超空化泡流动理论及应用[M].北京: 国防工业出版社, 2009.
[2]	李贝贝, 张宏超, 倪晓武, 等. 不同环境压强下激光空泡溃灭射流的实验研究[J].激光技术, 2012, 36(6):749–753.
[3]	马建旭, 王立鼎, 吴一辉. 微电子机械系统在生物医学领域中的应用[J].光学精密工程, 1996, 4(1):1–6.
[4]	杨拥军. 划时代的MEMS加工技术[J].国防制造技术, 2009, 10(5):36–42.
[5]	Mishra C, Peles Y. Cavitation in flow through a micro-orifice inside a silicon microchannel[J].Physics of Fluids, 2005, 17(1):151–155.
[6]	Schneider B, Ko?ar A, Kuo C J, et al. Cavitation enhanced heat transfer in microchannels[J].Journal of Heat Transfer, 2006, 128(12):1293–1301.DOI:10.1115/1.2349505
[7]	Mishra C, Peles Y. Flow Visualization of cavitating flows through a rectangular slot micro-orifice ingrained in a microchannel[J].Physics of Fluids, 2005, 17(11):3602–3614.
[8]	Mishra C, Peles Y. An experimental investigation of hydrodynamic cavitation in micro-venturis[J].Physics of Fluids, 2006, 18(10):3603–3605.
[9]	Fogg D W, Goodson K E. Bubble-induced water hammer and cavitation in microchannel flow boiling[J].Journal of Heat Transfer, 2009, 131(12):315–320.
[10]	Hickel S, Mihatsch M, Schmidt S J. Implicit large eddy simulation of cavitation in micro channel flows[J].Eprint Arxiv, 2014, 132(1):743–750.
[11]	Rooze J, André M, van der Gulik G S, et al. Hydrodynamic cavitation in micro channels with channel sizes of 100 and 750 micrometers[J].Microfluidics and Nanofluidics, 2012, 12(1-4):499–508.DOI:10.1007/s10404-011-0891-5
[12]	Schnerr G H, Saver J. Physical and numerical modeling of unsteady cativation dynamics[C]//Fourth International Conference on Multiphase Flow. New Orleans, Louisiana, USA, 2001.
[13]	Wilcox D C. Reassessment of the scale-determining equation for advanced turbulence models[J].AIAA Journal, 1988, 26(11):1299–1310.DOI:10.2514/3.10041
[14]	孙润鹏, 朱卫兵, 徐凌志, 等. 应用Transition k-kl-ω转捩模型对内冷叶片气热耦合的数值模拟[J].推进技术, 2012, 33(2):274–282.
[15]	Liu B, Cai J, Huai X L. Experimental study on heat transfer with cavitating flow in copper-based microchannels[C]// The Heat Transfer Symposium 2014. Beijing, China, 2014.