为保证原油管道的安全经济运行，须全面把握原油管道停输再启动过程，许多学者^[1-3]开展了胶凝原油管道测试研究。室内环道是常见测试手段之一，其测量原理是Hagen-Poiseuille流动。在Hagen-Poiseuille流动的初始阶段，必然存在管内流速分布不断变化的非稳态过程，造成测量黏度值与真实黏度值的差异。为明晰相关测试影响因素，减小测试误差，有必要掌握非稳态Poiseuille流动阶段流体表观的测试结果与真实流变性的差异以及不同边界条件对于非稳态流动过程的影响。最经典的非稳态Poiseuille流动为“置有牛顿流体的无限长水平圆管，初始静止，然后突然施加恒定压力，流体开始运动”。Papanastasiou等^[4]采用分离变量法给出了恒压力边界条件下圆管内不可压缩非稳态Poiseuille流动的解析解，但并未讨论恒流量条件下的非稳态Poiseuille流动。Barış^[5]对两种不可压缩牛顿流体混合物展开研究，得出非定常平板和轴对称Poiseuille流动下的速度分布解析解。Erdogan^[6-7]控制恒压力边界条件，对两个无限大平板间以及圆管内牛顿流体的非稳态Poiseuille流动展开研究。Erdogan等^[8]，研究了牛顿流体和二阶流体的剪切应力随时间变化规律的差异。Siddique等^[9]控制与时间有关的剪切应力边界，对不可压缩的广义Oldroyd-B流体进行研究，得到级数形式的速度分布和剪切应力的通用解。Fetecau^[10]对几何形状为轴对称的圆管的无边界区域内的Oldroyd-B流体进行研究。Muzychka等^[11]通过寻求非稳态Poiseuille流动和斯托克斯第一问题的关联，建立了描述内表面剪切应力以及速度分布随时间变化的简易模型，并对结果进行了渐近分析。Hayat等^[12-13]控制周期性压力边界和恒压力边界，对两平板间不可压缩二阶流体进行研究。此外，还有一些学者采用数值计算的思想展开了相关探究。Vinay等^[14]建立了非牛顿流体圆管内非稳态流动的数值计算模型，对方程组进行离散求解。Negrão等^[15]考虑了管道径向速度的变化开展研究。但上述学者仅讨论了恒压力边界条件下的非稳态流动过程。目前，关于Poiseuille流动的非稳态过程研究，主要关注流体非稳态速度分布函数的求解，边界条件集中在恒压力边界条件，并未关注不同边界条件下流动过程的差异以及非稳态过程对细管法测试黏度准确性的影响。笔者对不同边界条件下的圆管内非稳态Poiseuille流动过程进行分析，对比边界条件对于非稳态过程的影响。暂不考虑流体本身流变性的变化，选取牛顿流体为研究对象，考虑径向速度分布，采用数值计算方式，以无量纲黏度和无量纲时间为非稳态过程表征量，讨论不同边界条件对非稳态Poiseuille流动过程的影响。

1 问题描述 1.1 细管法黏度测试基础原理

细管法测定液体的流动性广泛应用于各种流体流变性的测量^[16]。对于稳定的管道层流，两端压差P恰好提供流体圆柱体在管道内壁上的黏滞力，管壁剪切应力τ与压差Δp之间满足

$ \tau = \frac{{\Delta pR}}{{2L}}. $

(1)

式中，L为细管长度，m；R为细管半径，m。

对于牛顿流体与幂律流体的稳定圆管层流，其管壁处剪切速率可表示为

$ \dot \gamma = \frac{Q}{{{\rm{ \mathsf{ π} }}{R^3}}}\frac{{3{n_{\rm{r}}} + 1}}{{{n_{\rm{r}}}}}. $

(2)

式中，$ \dot \gamma $为管壁剪切速率，s^-1；Q为管道流量，m³·s^-1；n_r为流变指数，当n_r=1时，为牛顿流体。

根据式(1)与(2)可获得流体剪切速率与剪切应力的关系，进而确定流体黏度或本构方程。测量黏度的表达式为

$ {\mu _{\rm{M}}} = \tau /\dot \gamma . $

(3)

式中，μ_M为测量黏度，Pa·s。

细管法开展流体测试，必须满足的核心条件是流体应为与时间无关的不可压缩流体，流动状态为充分发展的稳定层流^[16]。然而在实际流动过程中，在非稳态Poiseuille流动阶段，管内流体速度不断变化，基于细管法黏度测试原理计算的黏度数值不再是流体真实黏度。有必要明确非稳态Poiseuille流动阶段测试的流体表观流变性与理想均匀剪切条件下流体的真实流变性间的差异，从而探究初始非稳态流动阶段速度分布的发展对细管法测试的初始流变数据的影响。

1.2 表征非稳态过程的无量纲参数定义 1.2.1 无量纲黏度

在实际测试过程中的初始非稳态Poiseuille流动阶段，应用细管法准则计算的流体表观流变性与理想均匀剪切条件下流体的真实流变性存在差异，为表征这一差异，采用无量纲黏度$\tilde \mu $表征非稳态Poiseuille流动对于黏度测试结果的影响^[17-18]，$ \tilde \mu $的定义为

$ \tilde \mu = {\mu _{\rm{M}}}/\mu . $

(4)

式中，μ为牛顿流体真实黏度，Pa·s。当μ_M=μ时，$ \tilde \mu $=1，流动达到稳定，测量结果即为流体真实黏度。

使用无量纲黏度$ \tilde \mu $表征牛顿流体的非稳态Poiseuille流动特性，无量纲黏度越接近1，流动越趋于稳定。

1.2.2 无量纲时间

在利用解析方法进行非稳态流动的相关研究中通常利用无量纲时间描述非稳态流动过程^[7]，无量纲时间的定义为

$ \tilde t = \frac{t}{{\left( {{R^2}/\upsilon } \right)}} = \frac{{\upsilon t}}{{{R^2}}}. $

(5)

式中，$ \tilde t $为无量纲时间；υ为牛顿流体的运动黏度，m²·s^-1。

恒压力边界下非稳态过程中无量纲黏度$ \tilde \mu $仅是无量纲时间的函数，而在本文所选取的其他边界条件下非稳态过程的参数也可能仅与无量纲时间相关。将对不同边界条件下的非稳态过程进行无量纲化描述。

1.3 不同边界条件选取说明

在实际管道非稳态流动测试中恒流量^[19]和恒压力^[20-21]是最为常用的加载条件。除此之外，在控制流量边界下初期可能存在流量从0线性增加的情况^[22]。选取的边界条件分别为恒流量边界、流量从0线性增加边界和恒压力边界。

定义截面上平均速度概念为

$ \bar v\left( t \right) = \frac{Q}{{{\rm{ \mathsf{ π} }}{R^2}}} = \frac{{\int_0^R {2{\rm{ \mathsf{ π} }}rv{\rm{d}}r} }}{{{\rm{ \mathsf{ π} }}{R^2}}}. $

(6)

式中，$ \bar v\left( t \right) $为流体截面上的平均速度，m·s^-1；v为任意位置处微元流速，m·s^-1。

为计算方便所有的恒流量条件转换为恒平均速度边界条件，从0线性增加的流量条件转换为从0线性增加的平均速度边界。三大边界条件为:①恒平均速度边界条件；②从0线性增加的平均速度边界条件；③恒压力边界条件。

2 非稳态Poiseuille流动模型建立

对于不同边界条件下的圆管内Poiseuille流动，如图 1所示，取一段管长为L，管径为R的圆管，管道保持水平，管内充满不可压缩牛顿流体。初始状态全线压力、速度均为0。在管道始端突然施加一定的加载条件，流体开始运动。

如图 2所示，圆管内空间被划分为n个圆筒微元，i表示微元界面位置，I表示微元中心位置。设i=1等价于I=1，i=n等价于I=n。

2.1 动量守恒方程

对于径向第I(1 < I≤n)个圆筒微元，界面i与i-1合力作用的结果使其加速运动，动量方程为

$ \rho \frac{{v_I^{t + \Delta t} - v_I^t}}{{\Delta t}} = \frac{{\Delta p}}{L} + \left[ {\tau _{i - 1}^t\frac{{2\left( {i - 1} \right)}}{{2i - 1}} - \tau _i^t\frac{{2i}}{{2i - 1}}} \right]\frac{1}{{\Delta r}},1 \le I \le n. $

(7)

式中，上标t+Δt表示下一时刻，s；t表示当前时刻，s；ρ为流体密度，kg·m^-3；Δr为微元宽度，m；Δt为时间步长，s。

2.2 剪切应力

对于牛顿流体，其流变特性遵从牛顿内摩擦定律，即认为当流体受到外力作用时，其剪切速率与剪切应力的响应成正比。用流变方程描述为

$ \tau = \mu \dot \gamma . $

(8)

2.3 剪切速率

(1) I=0。界面处剪切速率可近似描述为

$ {{\dot \gamma }_{i = 0}} = 2\frac{{{v_{I = 0}} - {v_{I = 1}}}}{{\Delta r}} = 0. $

(9)

(2) 0 < I < n。界面处剪切速率可描述为

$ {{\dot \gamma }_i} = \frac{{{v_I} - {v_{I + 1}}}}{{\Delta r}}. $

(10)

(3) I=n。界面处剪切速率可描述为

$ {{\dot \gamma }_{i = n}} = \frac{{{v_{I = n}} - {v_{I = n + 1}}}}{{0.5\Delta r}} = 2\frac{{{v_{I = n}} - {v_{I = n + 1}}}}{{\Delta r}}. $

(11)

2.4 边界条件

(1) 恒平均速度边界。

$ \bar v\left( t \right) = {v_0}. $

(12)

式中，v₀为施加的恒定平均速度，m·s^-1。

(2) 从0线性增加的平均速度边界。

$ \bar v\left( t \right) = at,t > 0. $

(13)

式中，a为加速度，m·s^-2。

(3) 恒压力边界。

恒压力边界条件下圆管进出口压力均为恒定值，分别表示为

$ p\left( {x = 0} \right) = {p_0}, $

(14)

$ p\left( {x = L} \right) = 0. $

(15)

式中，x为圆管轴向坐标，m；p(x=0)表示圆管入口位置处的压力，Pa；p(x=L)为管道出口压力，Pa；p₀为施加的恒定压力，Pa。

2.5 初始条件

初始时刻，流体静止，则有

$ v_I^{t = 0} = 0,1 \le I \le n. $

(16)

剪切速率初始值为

$ {{\dot \gamma }^t} = {0_i} = 0,1 \le i \le n. $

(17)

剪切应力初始值为

$ {\tau ^t} = {0_i} = 0,1 \le i \le n. $

(18)

3 牛顿流体非稳态Poiseuille流动 3.1 数值计算基础参数定义

采用数值计算方法，通过编译程序对牛顿流体非稳态Poiseuille流动特性进行研究，采用控制单一变量法研究不同因素对于非稳态过程的影响规律。在数值计算过程中所需基础参数的默认值设置：μ=1 Pa·s，L=1 m，R=0.025 m，ρ=800 g·m^-3，Δt=10^-7 s，n=20。研究某一参数影响时，仅改变该参数数值，其他参数保持默认值不变，对比计算结果，分析该参数对于非稳态过程的影响规律。

3.2 数值计算模型验证

文献[4]中给出了恒压力边界条件下的牛顿流体非稳态Poiseuille流动速度分布函数的解析解。Erdogan ^[7]根据速度分布函数解析解，推导出流量随无量纲时间变化的解析解，表示为

$ \frac{{8QL\mu }}{{ - {\rm{ \mathsf{ π} }}\Delta p{R^4}}} = 1 - 32\sum\limits_{n = 1} {\frac{1}{{\lambda _n^4}}{{\rm{e}}^{ - \lambda _n^2\tilde t}}} . $

(19)

式中，Q为非稳态过程测量流量；-π(ΔP/Δr)R⁴/8μ为流动达到稳定后的流量；λ_n为贝塞尔函数。

圆管流动中，在恒压力边界条件下，黏度与流量成反比^[16]。即无量纲黏度应与无量纲流量互为倒数，则有

$ \tilde \mu = \frac{1}{{1 - 32\sum\limits_{n = 1} {\frac{1}{{\lambda _n^4}}{{\rm{e}}^{ - \lambda _n^2\tilde t}}} }}. $

(20)

通过编译程序，对牛顿流体进行数值求解，与解析解进行对比，验证数值计算模型的可靠性。无量纲黏度随无量纲时间的变化如图 3所示。

由图 3可知，由数值模型计算求解的无量纲黏度变化规律与文献[7]一致，说明所建立的模型具有一定的可靠性。

3.3 牛顿流体非稳态Poiseuille流动特点

根据上述牛顿流体非稳态Poiseuille流动的数值模型，编译程序，计算不同边界条件下的非稳态流动过程相关参数。

在Poiseuille流动的初始非稳态阶段，牛顿流体的速度分布随时间发生变化。以恒平均速度边界条件为例，取平均速度v₀为0.05 m·s^-1，其他参数根据默认值进行数值计算，记录不同时刻牛顿流体在圆管内的速度分布，如图 4所示。

从图 4可以看出，在流动发展的初始阶段，同一截面各位置处的流速相差较小，随着时间增加流速逐渐发展至稳定状态，管中心处流速最大，管壁处流速最小，整个截面中流速呈抛物线型分布，且不再随时间变化，仅与截面上所处位置有关，即流动得到充分发展。

由分析可知，Poiseuille流动的初始非稳态阶段主要是由径向速度分布发展的非稳态过程导致。

4 各因素对非稳态过程的影响规律 4.1 流体性质与管径参数影响 4.1.1 不同流体黏度下非稳态流动差异

计算牛顿流体非稳态Poiseuille流动过程中不同流体黏度条件下无量纲黏度$ \tilde \mu $随时间的变化规律，结果如图 5所示。

从图 5可以看出，随着时间增加无量纲黏度逐渐降至1，即测试黏度逐渐接近真实黏度。根据无量纲黏度曲线的变化规律可知，同一时刻的真实黏度越大，其对应的无量纲黏度越小，即测量偏差越小，无量纲黏度趋近1所需的时间也越短，即流动达到稳定所需的时间也越短。

4.1.2 不同管径下的非稳态流动差异

令管道直径D分别等于0.05、0.08、0.10 m，保持其他参数不变，分别计算牛顿流体非稳态Poiseuille流动过程无量纲黏度$ \tilde \mu $随时间变化规律。结果如图 6所示。

从图 6可以看出，随着时间增加无量纲黏度逐渐降至1，即测试黏度逐渐接近真实黏度。同一时刻对应的管径越大，其对应的无量纲黏度越大，即测量偏差越大，无量纲黏度趋近1所需的时间也越长，达到稳定流动所需的时间也越长。

4.1.3 非稳态阶段无量纲化描述

在恒压力边界下非稳态参数(如无量纲$ \tilde \mu $黏度)仅是无量纲时间的函数^[7]，而恒流量边界与线性增加流量边界条件下，非稳态过程的参数也可能仅与无量纲时间有关。

图 7为不同黏度和管径条件下无量纲黏度随无量纲时间变化。由图 7可知，不同黏度和管径条件下，牛顿流体在非稳态过程中的无量纲黏度随无量纲时间的变化规律完全一致，即无量纲黏度可视为仅与无量纲时间相关的函数。

4.2 不同边界条件下非稳态Poiseuille流动 4.2.1 边界条件参数取值的影响

分别在各边界条件下取不同数值进行计算。具体取值情况为：①恒平均速度边界条件下，取平均速度分别为0.01、0.05、0.10 m·s^-1；②从0线性增加的平均速度边界条件下，取加速度分别为0.01、0.1、1 m·s^-2；③恒压力边界条件下，别为600、1 000、2 000 Pa。

对比各边界条件，不同参数取值下无量纲黏度$ \tilde \mu $随无量纲时间$ \tilde t $的变化规律如图 8所示。

从图 8可以看出，不同速度、加速度及压力条件下，无量纲黏度$ \tilde \mu $随无量纲时间$ \tilde t $变化的曲线几乎重合。说明$ \tilde \mu $仅与$ \tilde t $有关，与平均速度边界的具体速度、从0线性增加的平均速度边界的具体加速度、恒压力边界的具体压力均无关。

4.2.2 边界条件类型影响

3种边界条件下的无量纲黏度对比如图 9所示。

从图 9可以看出，非稳态Poiseuille流动过程中$ \tilde \mu $随$ \tilde t $的增大而减小，并逐渐趋于1。其中，恒平均速度边界条件下达到$ \tilde \mu $=1所需的时间最短，恒压力边界条件次之，从0线性增加的平均速度边界条件所需时间最长。不同边界条件下Poiseuille流动达到均匀剪切所需的时间存在差异，为对差异进行进一步比较，给出定量描述：

$ \tilde \mu = \left( {\tilde t = {{\tilde t}_{\rm{d}}}} \right) = 1.01. $

(21)

式中，$ {\tilde t_{\rm{d}}} $为Poiseuille流动达到稳定对应的无量纲时间。

由于在非稳态流动阶段$ \tilde \mu $逐渐下降至1，所以定义当$ \tilde \mu $降至1.01，即测量黏度与真实黏度相对偏差等于1%时，流动达到稳定状态。根据数值计算模型计算得到3种边界条件Δp(t)=p₀、v(t)=v₀和v(t)=at(t>0)下的$ {\tilde t_{\rm{d}}} $分别为0.76、0.14和13.42。

相比恒压力边界，恒平均速度边界条件下Poiseuille流动达到稳定所需的时间更短。而相比恒平均速度边界，非定常的平均速度边界会造成Poiseuille流动的非稳态过程显著延长。

5 结论

(1) 3种边界条件下非稳态过程使得黏度测量结果较流体真实黏度偏大。随着非稳态阶段时间增大测量黏度逐渐趋近于真实黏度。

(2) 以无量纲黏度表征非稳态阶段特征，随着时间延长无量纲黏度逐渐降至1，且不同边界条件下达到1所对应的无量纲时间为定值。

(3) 在同一类型边界条件下无量纲黏度仅与无量纲时间有关，与边界条件中的平均速度、压力及加速度均无关。

(4) 对比不同类型边界条件，从0线性增加的平均速度边界会造成Poiseuille流动的非稳态过程显著延长，恒压力边界条件所需的流动发展时间次之，恒平均速度边界条件对应的非稳态过程最短。