清华大学 热能工程系, 热科学与动力工程教育部重点实验室, 北京 100084
收稿日期：2016-03-23
作者简介：宁涛(1989—), 男, 博士研究生
通信作者：顾春伟, 教授, E-mail:gcw@tsinghua.edu.cn
摘要：该文利用商业软件ANSYS-CFX对圆管中的脉动湍流进行了大涡模拟，分析了脉动流的减阻特性和总能耗。文中的脉动流算例包括稳态流主控和振荡流主控两种流态。结果表明：脉动流通过叠加适合的振荡流来改变稳态流的边界层特性，脉动幅值为5.5时得到最佳减阻率为25%；当脉动流的流态由振荡流主控且振荡流分量的边界层为层流时，减阻效果较好；简单正弦形式脉动流的总能耗高于相应的稳态流。
关键词：圆管脉动流湍流大涡模拟减阻能耗
Numerical analysis of the drag reduction for turbulent pulsating pipe flows based on large eddy simulations
NING Tao, GU Chunwei
Key Laboratory for Thermal Science and Power Engineering of Ministry of Education, Department of Thermal Engineering, Tsinghua University, Beijing 100084, China
Abstract: Large eddy simulations were conducted for turbulent pulsating flows using the commercial solver ANSYS-CFX. The drag reduction and the total energy consumption for pulsating flows were analyzed. The simulations included current dominated and wave dominated pulsating flows. The boundary layer characteristics of the current flow were affected by the superposition of the wave flow. The best drag reduction in the pulsating flows gave a 25% drag reduction when the non-dimensional pulsating amplitude was 5.5. The analysis indicates that the drag reduction is optimized when the pulsating flow is wave dominated and the wave boundary layer is laminar. Pulsating flows with simple sinusoidal pulsating patterns consume much more energy than steady flows.
Key words: pulsating pipe flowsturbulencelarge eddy simulationdrag reductionenergy consumption
脉动流是一种典型的非稳态流动，可以看作由稳态流和时均速度为零的振荡流叠加而成，其流动参数如压力、速度等会产生周期性变化。相比于稳态流，脉动流增加了一个时间尺度(振动频率)和一个空间尺度(振动幅值)，使得对脉动流的研究异常艰难[1]。到目前为止，脉动流的流动机理尚不清楚，已有研究仅限于圆管和槽道中的简单正弦脉动流。脉动流的流动特征可以用稳态Reynolds数(current Reynolds number, ReC)、振荡Reynolds数(wave Reynolds number, ReW)和脉动幅值(α)这3个量纲为1的参数来描述[1]。以圆管为例，其定义如下：
| ${{\mathop{\rm \mathit{Re}}\nolimits} _{\rm{C}}} = VD/{U_{\rm{m}}}, {{\mathop{\rm \mathit{Re}}\nolimits} _{\rm{W}}} = U_{\rm{m}}^2/\left( {\omega v} \right), \alpha = {U_{\rm{m}}}/V.$ | (1) |
本文以圆管脉动湍流为研究对象，采用大涡模拟(large eddy simulation, LES)方法，选用Dynamic Smagorinsky亚格子应力模式，主要对脉动湍流的减阻特性进行了研究，并对脉动湍流的总耗能进行了评估, 探讨脉动湍流在管道输运等实际工程中的应用潜力。
1 研究对象本文以水平光滑圆管为研究对象，以不可压缩的水为流动介质。圆管直径为D=2R=0.2 m，圆管长度为L=4πR。
由式(1) 可知，圆管脉动流的参数可以由ReC、ReW和α这3个量纲为1的参数来描述。本文共设置了5组算例，脉动频率为固定值，表征频率的Stokes数取为R/δ=53，各算例的其他参数如表 1所示。
表 1 稳态流和脉动流算例的相关参数
| 算例 | ReC | ReW | α |
| 稳态流 | 6 087 | 0 | 0 |
| 脉动流1 | 6 077 | 103 | 0.25 |
| 脉动流2 | 6 116 | 2 470 | 1.3 |
| 脉动流3 | 6 834 | 62 629 | 5.5 |
| 脉动流4 | 6 019 | 280 321 | 13 |
表选项
2 数值模拟方法2.1 数值方法Scotti等[8]采用大涡数值模拟对槽道脉动湍流的研究中，通过与直接数值模拟(direct numerical simulation, DNS)以及实验的对比，确认了Dynamic Smagorinsky亚格子应力模式能够很好地预测脉动湍流。故本文采用ANSYS软件包中的CFX作为求解器，选用Dynamic Smagorinsky亚格子应力模式对管流进行大涡数值模拟。其中:扩散项离散采用中心型格式，时间项离散采用二阶Euler后差格式。
2.2 边界条件本文假设圆管周向和流向的湍流统计存在周期性，为得到充分发展的湍流，将圆管进口和出口设为位移周期性边界条件。针对稳态流算例，进出口给定固定压强梯度Δp/L；针对脉动流算例，进出口给定随时间正弦变化的压强梯度(Δp/L)·1+βsin(ωt)。其中: β与脉动幅值相关，ω表征脉动频率。本文中Stokes数取为R/δ=53，对应的频率取为固定值ω=0.501，压差Δp=0.08 Pa。
2.3 网格划分本文对圆管计算域中的网格拓扑划分采用里层H型网格、外层O形网格的形式，如图 1所示。计算网格总数为920 623，流向、径向、周向的量纲为1的网格尺度如表 2所示，其定义如下：
| $\left\{ {\begin{array}{*{20}{c}}{{\rm{流向}}\;\;\;\;\;\;\;\;\;\;\;\;\;\;\Delta {z^ + } = \Delta z{u_\tau }/v, }\\{{\rm{径向}}\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\Delta {r^ + } = \Delta r{u_\tau }/v, }\\{{\rm{周向}}\;\;\;\;\;\;\;{{\left( {r\Delta \theta } \right)}^ + } = \left( {r\Delta \theta } \right){u_\tau }/v.}\end{array}} \right.$ | (2) |
 |
| 图 1 计算网格示意 |
| 图选项 |
表 2 量纲为1的网格尺度和时间尺度
| 网格总数 | Δz+ | Δr+ | (rΔθ)+ | Δt+ |
| 920 623 | 20.0 | 0.02~6.6 | 6.6~12.6 | 7×10－5 |
表选项
其中：
计算物理时间步长取为Δt=0.008 s，对应的量纲为1的时间步长为Δt+=Δtuτ/D=7×10－5。
3 计算结果及分析3.1 计算验证本文假设湍流统计量沿流向和周向存在周期性，将各流动参数先沿流向和周向做空间平均，再进一步做时间平均，得到流动参数的系综统计平均值。
为了验证求解器CFX中的LES对本文研究问题预测的准确性，本文将稳态流算例的数值计算结果与文[11-12]中实验研究结果以及文[12]中的直接数值模拟结果进行对比。
稳态流计算的Reynolds数为6 087。下面分别对沿径向的量纲为1的截面平均速度U+和湍动度T进行对比，U+和T定义如式(3) 所示：
| ${U^ + } = V/{u_\tau }, T = \sqrt {u'u'} /{u_\tau }.$ | (3) |
图 2对比了本文稳态流算例和文[12]中沿径向的量纲为1的速度剖面。在y+ < 5的黏性底层区域，本文计算结果是速度分布满足U+=y+的线性规律；在y+>12的对数层区域，速度分布稍微偏离对数分布U+=2.5 lny++5.5，但是与文[12]中的实验值符合良好。
 |
| 图 2 量纲为1的速度剖面 |
| 图选项 |
图 3给出了沿径向的湍动度分布，其分布趋势与文[11-12]中的实验值吻合良好。与文[12]中的直接数值模拟相比，湍动度T均在y+=14附近达到U+=y+极大值。由于Reynolds数略有差别，直接数值模拟的湍动度极大值略微偏小。
 |
| 图 3 沿径向的湍动度对比 |
| 图选项 |
通过对平均速度剖面和湍动度的对比，CFX中的Dynamic Smagorinsky亚格子应力模式能够比较精确地预测本文的研究问题。
3.2 减阻分析为研究脉动湍流的减阻特性，本文共计算了4个脉动湍流的算例，减阻率DR定义为
| ${D_{\rm{R}}} = \left( {{C_{{\rm{fBlasius}}}}-{C_{\rm{f}}}} \right)/{C_{{\rm{fBlasius}}}}.$ | (4) |
减阻效果如表 3所示。可见，脉动幅值α=5.5对应算例的减阻效果最好，减阻率DR达到25%，下面对该算例进行详细分析。
表 3 不同脉动流算例的参数和减阻率
| 算例 | ReC | ReW | α | DR/% |
| 脉动流1 | 6 077 | 103 | 0.25 | 5 |
| 脉动流2 | 6 116 | 2 470 | 1.3 | 10 |
| 脉动流3 | 6 834 | 62 929 | 5.5 | 25 |
| 脉动流4 | 6 091 | 280 321 | 13 | 18 |
表选项
图 4给出了脉动幅值α=5.5对应算例在φ分别为0、π/2、π、3π/2这4个不同相位的相平均速度剖面。相平均即各个相位所对应的系综平均。由于脉动幅值α很大，可以看出不同相位的相平均速度剖面变化较大，但是由于Stokes数R/δ=53很大，不同相位的对数层剖面形状基本不变。
 |
| 图 4 不同相位的相平均速度剖面(α=5.5) |
| 图选项 |
图 5给出了α=5.5对应的脉动流3算例和稳态流算例的Reynolds应力
 |
| 图 5 稳态流和脉动流的Reynolds应力 |
| 图选项 |
图 6给出了4个脉动流算例的周期平均的速度剖面。可以看出，4个算例的对数层区域的斜率基本不变，相比于稳态流算例，脉动流算例的对数层有不同程度的上移。对数层的上移间接表明黏性底层厚度的增加，而黏性底层厚度的增加是流动减阻的标志性特征。α=5.5对应的脉动流的对数层上移最多，对应于最大的减阻效率，而其他算例的对数层上移程度和减阻率的大小也是互相对应的。
 |
| 图 6 不同算例的速度剖面比较 |
| 图选项 |
下面进一步讨论脉动流减阻的最佳参数区间。根据脉动幅值α可以将脉动流的流态分为稳态流主控和振荡流主控两种，根据振荡流Reynolds数ReW的临界值1.5×105可以将振荡流的边界层分为层流和湍流[1]: ReW < 1.5×105时，振荡流边界层为层流；ReW>1.5×105时，振荡流边界层为湍流。图 7给出了本文和文[9]中大涡数值模拟得到的减阻效果。可以看出，当脉动流中叠加的振荡流分量的边界层为层流时，随着脉动幅值的增大，减阻率增加; 当振荡流分量的边界层为湍流时，随着脉动幅值的增大，减阻率降低。脉动流是通过叠加合适的振荡流来改变稳态流边界层的特性从而达到减阻的目的。当叠加的振荡流的边界层为层流且振荡流分量主控脉动流流态时，稳态湍流分量有被振荡流分量层流化的趋势。从计算结果和以上分析可以得到脉动流减阻的最佳参数区间：脉动流的流态由振荡流主控，且振荡流分量的边界层为层流。
 |
| 图 7 大涡模拟得到的减阻率比较 |
| 图选项 |
3.3 总能耗评估稳态流叠加适当的振荡流可以减小流动阻力，但是产生振荡流需要消耗额外的能量。先前的研究[1-12]一般只针对脉动流的减阻特性，很少有研究者评价脉动流的总能耗。对实际工程应用而言，脉动流的总能耗评估是必要的，总能耗可以用计算周期平均的做功量来表示，本文采用式(5) 对稳态流和脉动流的能量消耗进行评价：
| $\left\{ {\begin{array}{*{20}{c}}{{E_{\rm{P}}} = \frac{{\int_0^{2{\rm{\pi }}} {u\left( \phi \right)\Delta p\left( {1 + \beta \sin \phi } \right)R{\rm{d}}\phi } }}{{2{\rm{\pi }}\rho L{V^3}}}, }\\{{E_{\rm{S}}} = \frac{{\int_0^{2{\rm{\pi }}} {V\Delta pR{\rm{d}}\phi } }}{{2{\rm{\pi }}\rho L{V^3}}} = \frac{{\Delta pR{\rm{d}}\phi }}{{\rho L{V^2}}}.\;\;\;\;\;\;}\end{array}} \right.$ | (5) |
不同脉动流算例的总能耗如表 4所示。可以看出，脉动流和相应稳态流的能耗比EP/ES都大于1，表明脉动流的总能耗高于相应稳态流，且总能耗随着脉动幅值α的增大而增大。这说明本文研究的正弦形式脉动流虽然在一定参数范围内能使平均阻力显著降低，但是总能耗是增加的。最近Souma等[13]在实验研究中发现在脉动周期中改变加速阶段和减速阶段的频率和幅值可以使总能耗减小。可见，脉动减阻作为实际工程应用中的一种减阻手段还需要进一步研究。
表 4 不同算例的脉动总能耗比较
| 算例 | α | EP | ES | EP/ES |
| 脉动流1 | 0.25 | 0.009 8 | 0.008 | 1.13 |
| 脉动流2 | 1.3 | 0.035 6 | 0.008 | 4.15 |
| 脉动流3 | 5.5 | 0.472 7 | 0.006 | 68.90 |
| 脉动流4 | 13.0 | 3.087 4 | 0.008 9 | 349.0 |
表选项
4 结论本文利用商业软件ANSYS-CFX对圆管中的脉动湍流进行了大涡模拟，主要研究了脉动湍流的减阻特性，并对其总耗能进行了评估，得到以下结论:
1) 本文计算的脉动流算例包括稳态流主控和振荡流主控两种流态，在脉动幅值α=5.5时获得了显著的减阻效果，减阻率DR为25%。
2) 脉动流主要通过叠加合适的振荡流来改变稳态流边界层的特性从而达到减阻目的。当振荡流分量的边界层为层流且振荡流分量主控脉动流流态时，稳态湍流分量有被振荡流分量层流化的趋势。因此，当脉动流的流态由振荡流主控且振荡流分量的边界层为层流时，减阻效果较好。
3) 脉动流的总能耗高于相应的稳态流，采用非常规的脉动形式有减小总能耗的潜力。要将脉动流作为实际工程应用中的一种减阻手段还需要进一步研究。
参考文献
| [1] | Journal of Central South University(Science and Technology), 41(2):649-654.-->Lodahl C R, Sumer B M, Fredsoe J. Turbulent combined oscillatory flow and current in a pipe[J]. Journal of Fluid Mechanics, 1998, 373: 313–348. DOI:10.1017/S0022112098002559 |
| [2] | Journal of Central South University(Science and Technology), 41(2):649-654.-->MAO Zhuoxiong, Hanratty T J. Studies of the wall shear stress in a turbulent pulsating pipe flow[J]. Journal of Fluid Mechanics, 1986, 170: 545–564. DOI:10.1017/S0022112086001015 |
| [3] | Journal of Central South University(Science and Technology), 41(2):649-654.-->Akhavan R, Kamm R D, Shapiro A H. An investigation of transition to turbulence in bounded oscillatory stokes flows, Part 1. Experiments[J]. Journal of Fluid Mechanics, 1991, 225: 395–422. DOI:10.1017/S0022112091002100 |
| [4] | Journal of Central South University(Science and Technology), 41(2):649-654.-->Tu S W, Ramaprian B R. Fully developed periodic turbulent pipe flow: Part 1. Main experimental results and comparison with predictions[J]. Journal of Fluid Mechanics, 1983, 137: 31–58. DOI:10.1017/S0022112083002281 |
| [5] | Journal of Central South University(Science and Technology), 41(2):649-654.-->Ramaprian B R, Tu S W. Fully developed periodic turbulent pipe flow: Part 2. The detailed structure of the flow[J]. Journal of Fluid Mechanics, 1983, 137: 59–81. DOI:10.1017/S0022112083002293 |
| [6] | Journal of Central South University(Science and Technology), 41(2):649-654.-->He S, Jackson J D. An experimental study of pulsating turbulent flow in a pipe[J]. European Journal of Mechanics-B: Fluids, 2009, 28(2): 309–320. DOI:10.1016/j.euromechflu.2008.05.004 |
| [7] | Journal of Central South University(Science and Technology), 41(2):649-654.-->Tuzi R, Blondeaux P. Intermittent turbulence in a pulsating pipe flow[J]. Journal of Fluid Mechanics, 2008, 599: 51–79. |
| [8] | Journal of Central South University(Science and Technology), 41(2):649-654.-->Scotti A, Piomelli U. Numerical simulation of pulsating turbulent channel flow[J]. Physics of Fluids, 2001, 13(5): 1367–1384. DOI:10.1063/1.1359766 |
| [9] | Journal of Central South University(Science and Technology), 41(2):649-654.-->Manna M, Vacca A. Resistance reduction in pulsating turbulent pipe flows[J]. Journal of Engineering for Gas Turbines and Power, 2005, 127: 410–417. DOI:10.1115/1.1789511 |
| [10] | Journal of Central South University(Science and Technology), 41(2):649-654.-->Manna M, Vacca A. Spectral dynamic of pulsating turbulent pipe flow[J]. Computers & Fluids, 2008, 37(7): 825–835. |
| [11] | Journal of Central South University(Science and Technology), 41(2):649-654.-->Tardu S F, Binder G, Blackwelder R F. Turbulent channel flow with large-amplitude velocity oscillations[J]. Journal of Fluid Mechanics, 1994, 267: 109–151. DOI:10.1017/S0022112094001138 |
| [12] | Journal of Central South University(Science and Technology), 41(2):649-654.-->Eggels J G M, Unger F, Weiss M H, et al. Fully developed turbulent pipe flow: A comparison between direct numerical simulation and experiment[J]. Journal of Fluid Mechanics, 1994, 268: 175–210. DOI:10.1017/S002211209400131X |
| [13] | Journal of Central South University(Science and Technology), 41(2):649-654.-->Souma A, Iwamoto K, Murata A. Experimental analysis of pressure-gradient profile upon drag-reduction effect in pulsating turbulent pipe flow[J]. Transactions of the Japan Society of Mechanical Engineers Series B, 2012, 78(787): 521–530. DOI:10.1299/kikaib.78.521 |
