脉动流是一种典型的非稳态流动，可以看作由稳态流和时均速度为零的振荡流叠加而成，其流动参数如压力、速度等会产生周期性变化。相比于稳态流，脉动流增加了一个时间尺度(振动频率)和一个空间尺度(振动幅值)，使得对脉动流的研究异常艰难^[1]。到目前为止，脉动流的流动机理尚不清楚，已有研究仅限于圆管和槽道中的简单正弦脉动流。脉动流的流动特征可以用稳态Reynolds数(current Reynolds number, Re_C)、振荡Reynolds数(wave Reynolds number, Re_W)和脉动幅值(α)这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)

其中：D为圆管直径，V为时间统计平均的截面平均速度, U_m为圆管轴线上的速度脉动幅值，ω=2π/T为脉动角速度，ν为运动黏性系数。当脉动幅值α$ \ll $1时，脉动流的流态为稳态流主导；当α$\gg $1时，脉动流的流态为振荡流主导^[1]。迄今为止，关于圆管和槽道湍流脉动的研究大多仅局限于稳态流主导的脉动流^[2-8]，仅有少数研究涉及振动流主导的湍流脉动，例如Lodahl等^[1]的实验研究和Manna等^[9-10]的数值研究表明，在振荡流主导的一些参数区间内，周期平均的壁面切应力与稳态流相比有显著降低。但是，与稳态流相比，产生脉动流需要消耗额外的能量，实际工程应用(如石油、天然气管道输运的减阻节能)必须考虑总耗能。

本文以圆管脉动湍流为研究对象，采用大涡模拟(large eddy simulation, LES)方法，选用Dynamic Smagorinsky亚格子应力模式，主要对脉动湍流的减阻特性进行了研究，并对脉动湍流的总耗能进行了评估, 探讨脉动湍流在管道输运等实际工程中的应用潜力。

1 研究对象

本文以水平光滑圆管为研究对象，以不可压缩的水为流动介质。圆管直径为D=2R=0.2 m，圆管长度为L=4πR。

由式(1) 可知，圆管脉动流的参数可以由Re_C、Re_W和α这3个量纲为1的参数来描述。本文共设置了5组算例，脉动频率为固定值，表征频率的Stokes数取为R/δ=53，各算例的其他参数如表 1所示。

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)

其中：${u_\tau } = \sqrt {R\Delta p/\left( {2\rho L} \right)} $为壁面摩擦速度; 以z、r和θ分别表示圆管的流向、径向和周向坐标; 以u、v和w分别表示速度的流向、径向和周向分量。

计算物理时间步长取为Δ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)

其中：V表示系综平均的截面平均速度，u_τ表示时间和空间平均的壁面摩擦速度，u′表示沿流向的速度脉动。

图 2对比了本文稳态流算例和文[12]中沿径向的量纲为1的速度剖面。在y⁺ < 5的黏性底层区域，本文计算结果是速度分布满足U⁺=y⁺的线性规律；在y⁺>12的对数层区域，速度分布稍微偏离对数分布U⁺=2.5 lny⁺+5.5，但是与文[12]中的实验值符合良好。

图 3给出了沿径向的湍动度分布，其分布趋势与文[11-12]中的实验值吻合良好。与文[12]中的直接数值模拟相比，湍动度T均在y⁺=14附近达到U⁺=y⁺极大值。由于Reynolds数略有差别，直接数值模拟的湍动度极大值略微偏小。

通过对平均速度剖面和湍动度的对比，CFX中的Dynamic Smagorinsky亚格子应力模式能够比较精确地预测本文的研究问题。

3.2 减阻分析

为研究脉动湍流的减阻特性，本文共计算了4个脉动湍流的算例，减阻率D_R定义为

$ {D_{\rm{R}}} = \left( {{C_{{\rm{fBlasius}}}}-{C_{\rm{f}}}} \right)/{C_{{\rm{fBlasius}}}}. $

(4)

式中：C_f=2τ_W/(ρV²)表示计算得到的壁面摩擦系数，其中$ {\tau _{\rm{W}}} = \frac{1}{T}{\int_0^T {\tilde \tau } _{\rm{W}}}\left( {\omega t} \right){\rm{d}}t$表示系综周期平均的壁面切应力；C_fBlasius=0.079 Re_C^－0.25表示以脉动湍流的稳定流Reynolds数Re_C为基准的Blasius壁面摩擦系数。

减阻效果如表 3所示。可见，脉动幅值α=5.5对应算例的减阻效果最好，减阻率D_R达到25%，下面对该算例进行详细分析。

图 4给出了脉动幅值α=5.5对应算例在φ分别为0、π/2、π、3π/2这4个不同相位的相平均速度剖面。相平均即各个相位所对应的系综平均。由于脉动幅值α很大，可以看出不同相位的相平均速度剖面变化较大，但是由于Stokes数R/δ=53很大，不同相位的对数层剖面形状基本不变。

图 5给出了α=5.5对应的脉动流3算例和稳态流算例的Reynolds应力$\overline {u'v'} $沿径向的分布。稳态流和脉动流的Reynolds应力极大值都出现在y⁺=50附近，脉动流的Reynolds应力相比于稳态流沿整个径向都有较大程度的减小，这和该脉动流算例显著的减阻效果是一致的。

图 6给出了4个脉动流算例的周期平均的速度剖面。可以看出，4个算例的对数层区域的斜率基本不变，相比于稳态流算例，脉动流算例的对数层有不同程度的上移。对数层的上移间接表明黏性底层厚度的增加，而黏性底层厚度的增加是流动减阻的标志性特征。α=5.5对应的脉动流的对数层上移最多，对应于最大的减阻效率，而其他算例的对数层上移程度和减阻率的大小也是互相对应的。

下面进一步讨论脉动流减阻的最佳参数区间。根据脉动幅值α可以将脉动流的流态分为稳态流主控和振荡流主控两种，根据振荡流Reynolds数Re_W的临界值1.5×10⁵可以将振荡流的边界层分为层流和湍流^[1]: Re_W < 1.5×10⁵时，振荡流边界层为层流；Re_W>1.5×10⁵时，振荡流边界层为湍流。图 7给出了本文和文[9]中大涡数值模拟得到的减阻效果。可以看出，当脉动流中叠加的振荡流分量的边界层为层流时，随着脉动幅值的增大，减阻率增加; 当振荡流分量的边界层为湍流时，随着脉动幅值的增大，减阻率降低。脉动流是通过叠加合适的振荡流来改变稳态流边界层的特性从而达到减阻的目的。当叠加的振荡流的边界层为层流且振荡流分量主控脉动流流态时，稳态湍流分量有被振荡流分量层流化的趋势。从计算结果和以上分析可以得到脉动流减阻的最佳参数区间：脉动流的流态由振荡流主控，且振荡流分量的边界层为层流。

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)

式中：E_P表示脉动流的量纲为1的总能耗，E_S表示相应稳态流的量纲为1的总能耗，u(φ)表示相应φ相位的系综平均的截面平均速度。

不同脉动流算例的总能耗如表 4所示。可以看出，脉动流和相应稳态流的能耗比E_P/E_S都大于1，表明脉动流的总能耗高于相应稳态流，且总能耗随着脉动幅值α的增大而增大。这说明本文研究的正弦形式脉动流虽然在一定参数范围内能使平均阻力显著降低，但是总能耗是增加的。最近Souma等^[13]在实验研究中发现在脉动周期中改变加速阶段和减速阶段的频率和幅值可以使总能耗减小。可见，脉动减阻作为实际工程应用中的一种减阻手段还需要进一步研究。

4 结论

本文利用商业软件ANSYS-CFX对圆管中的脉动湍流进行了大涡模拟，主要研究了脉动湍流的减阻特性，并对其总耗能进行了评估，得到以下结论:

1) 本文计算的脉动流算例包括稳态流主控和振荡流主控两种流态，在脉动幅值α=5.5时获得了显著的减阻效果，减阻率D_R为25%。

2) 脉动流主要通过叠加合适的振荡流来改变稳态流边界层的特性从而达到减阻目的。当振荡流分量的边界层为层流且振荡流分量主控脉动流流态时，稳态湍流分量有被振荡流分量层流化的趋势。因此，当脉动流的流态由振荡流主控且振荡流分量的边界层为层流时，减阻效果较好。

3) 脉动流的总能耗高于相应的稳态流，采用非常规的脉动形式有减小总能耗的潜力。要将脉动流作为实际工程应用中的一种减阻手段还需要进一步研究。