水平同轴转柱体间隙区域间水动力数值模拟研究

谈雯倩，张东来，邵宇阳*，郝林杰，诸裕良

(1.河海大学 港口海岸与近海工程学院，南京 210098；2.中交上海航道勘察设计研究院有限公司，上海 201208；3.中水北方勘测设计研究有限责任公司，天津 300222)

细颗粒泥沙为主的河口海岸的航道往往会有泥沙淤积碍航的困扰，这一问题会严重影响相关流域的生产和经济发展。其中细颗粒泥沙的絮凝过程直接关系到河口地区的泥沙沉降-淤积状况[1]，是河口地形演变、航槽清淤治理、维护水体质量等问题中的一个重要的研究对象。前人针对絮凝过程开展了大量研究，其中絮凝实验常用的装置有搅动装置[2-3]、泰勒库特流装置[4-5]等，泰勒库特流装置因其操作便捷，通过改变同轴转柱间隙内流体的流速梯度来研究对泥沙絮凝的影响，控制精准，而被广泛使用。

流体在旋转的两个同轴圆柱之间的流动，即泰勒库特流，Couette[6]最早在流体粘性系数研究实验中对此展开研究，之后由Taylor[7]进一步研究该流动的稳定性理论，泰勒库特流是流体力学中一个重要的研究课题。国内外学者关于泰勒库特流的研究有很多，早期的研究主要是实验，Cloe[8]研究了不同半径比下的泰勒涡和波动泰勒涡出现所需的临界转速，且得出了产生波动泰勒涡的临界转速与内外圆柱的间隙大小成反比的结论。Werely[9]等人试验研究了仅内圆柱转动的情况下间隙内泰勒库特流。Akonur[10]等利用PIV粒子图像测速法研究了仅内圆柱旋转的条件下泰勒库特流的子午面速度，实验结果发现泰勒涡的强度随着雷诺数的增加而增强。鲍锋[11]等人通过流动显示实验和PIV粒子图像测速法对间隙流场进行可视化和定量化的研究，分析流体运动的周期性规律以及间隙内的雷诺数分布。近年来数值模拟逐渐成为常用的研究方法，Hwang[12]等人以Werely[9]的实验结果为研究基础，建立数学模型比较了环形间隙内有无轴向流动对流态的影响，验证了轴向流动对流场的稳定作用。冯俊杰[13]将实验和数值模拟相结合得出内圆柱转速、旋转雷诺数与流场涡形态转变的关系。Berghout[14]将Monin-Obukhov曲率长度的概念应用于泰勒库特流的研究，得到了湍流状态下的平面速度剖面。Cheng[15]等人利用大涡模拟捕捉泰勒库特流流动，并提出粗糙度模型以作为今后泰勒库特流设计小尺度均匀表面粗糙度实验的指导。泰勒库特流的数值模拟是近年来研究的热点问题，泰勒库特流模型在军工、航空、动力工程、生物工程、水处理、防洪、生态环境保护和膜分离等领域有着广泛的应用，因此对它的研究具有重要的理论和实践意义。

综合目前的研究来看，对泰勒库特流的研究一般为竖直放置，泥沙在竖直放置的泰勒库特流装置中受到重力的作用，会在轴向上存在一个泥沙浓度梯度，有研究证明泥沙浓度亦是絮凝研究的主要影响因素之一[16]。将同轴转柱水平放置，周向的流动可以使纵向上泥沙搅动均匀，即可有效避免此类影响。目前对水平的泰勒库特流研究较少，因此本文利用Fluent软件对三维不可压缩的、仅内圆柱旋转的泰勒库特流流场进行模拟研究，对比前人的实验结果进行模型验证，论证数值模拟的可行性，并在此基础上建立模型模拟水平放置的同轴转柱体间隙区域间流场，研究其水动力特征，为今后的絮凝实验研究打好基础。

1 模型建立

Wereley和Lueptow[17]曾经对泰勒库特流进行了实验研究，所使用的装置轴向长度为41 cm，内圆柱直径为8.68 cm，外圆柱直径为10.46 cm，间隙宽度为0.89 cm。装置内的液体为水、甘油和碘化钠的混合液，液体密度为1 620 kg/m3，运动粘性系数为3.15×10-6m2/s。在实验中，仅内圆柱旋转，通过改变内圆柱转速，研究泰勒库特流装置内的流场变化，得到了间隙流态由层流转变为层流泰勒涡的临界雷诺数与雷诺数为124时的层流泰勒涡流场。本文根据其实验结果建立模型，并进行模型验证。

模拟泰勒库特流装置，选择初始条件和边界条件，所建立的泰勒库特流模型流场由静止开始启动，内壁面为旋转壁面，以设定的转速均匀转动，外壁面为静止壁面。

1.1 网格模型

通过数值模拟计算控制方程确定解时，需要将控制方程离散化。Fluent中采用有限体积法对控制方程组进行空间和时间上的离散，然后使用分离式解法对离散后的方程组进行求解。

本文模拟的泰勒库特流装置，计算区域为环形，选择均匀分布的结构化六面体核心网格进行划分。模拟同轴圆柱之间的流体流动，需同时考虑径向、轴向和周向三个方向上的网格数量，本文在三个方向上划分的网格数量依次为32、128和256个。

在划分网格时，考虑到边界层的影响，在内、外壁面和上、下端面共四个壁面内设置了膨胀层，近壁面起始网格高度设置为0.1 mm，膨胀系数为1.1。

1.2 湍流模型

Realizablek-ε两方程模式适用于旋转剪切流、含有射流、混合流的自由流动、管道内流动、边界层流动、带有分离的流动和有明显旋流的流动。有对比分析认为Realizablek-ε湍流模型可以很好地模拟泰勒库特流的湍流状态[18]。Realizablek-ε两方程模式的湍流能量输运方程和能量耗散方程如下

(1)

(2)

式中：右端三项分别代表生成项、耗散项和壁面项。

2 模拟结果验证

Wereley等人的实验重点关注了泰特库特流在雷诺数为103以及124时的流动状态，本文将参照其实验结果，进行数值模拟计算。

2.1 临界雷诺数的确定

同轴圆柱间隙内流体由层流转变为层流泰勒涡时的雷诺数，称为临界雷诺数。Wereley等人的实验中，当内圆柱的转速达到0.827 rad/s时，流态发生变化，并得知此时的雷诺数为103。根据此数据，本文模拟雷诺数为103时的流场，发现因端部效应仅两端出现涡状，其余中间部分仍处于层流状态，总体流态基本可以认为是层流。加大转速，模拟雷诺数为104时的流动状态，开始出现了层流泰勒涡的特征，有微弱的涡状，而当雷诺数Re=105时，层流泰勒涡的特征更为明显，流场出现规则分布的涡状，间隙内的流动状态完全由层流转换为完全的层流泰勒涡。

2.2 临界雷诺数时流场模拟结果

2.2.1 临界雷诺数时流态分布

临界雷诺数子午面速度矢量图如图1所示。上实线代表泰勒库特流装置内壁，下实线为外壁，内部曲线为速度等值线，等值线之间的差值为0.1倍内壁转速。

图1 临界雷诺数子午面速度矢量图Fig.1 Velocity vector of meridian plane at critical Reynolds number

图1-a为Wereley等通过实验得到的结果，本文模拟所得的临界雷诺数在104和105附近，则此时的子午面流速矢量图如图1-b、1-c所示。

通过图1-b、图1-c可以看出，雷诺数为104时，水流泰勒涡状不明显，雷诺数为105时，泰勒涡充分发展，流态涡型明显。

子午面速度矢量图体现了径向流动造成的速度变化。图1-a、图1-c等值线变形较为一致，内壁处高速流体向外流动，外壁处低速流体向内流动，从而形成涡状。明显的流体向外流动的速度大于向内流动，因而向外流动区域的速度等值线的变化大于向内流动区域速度等值线的变化。

2.2.2 临界雷诺数时径向流动特性

图2为子午面间隙中轴线上的径向速度对比图。实验利用PIV粒子图像测速法测得结果。数值模拟雷诺数为104和105时的径向速度，其值与实验结果有一定的差距，进而模拟雷诺数为106时的结果，与实验结果相似。

图2 子午面中轴线径向速度对比图Fig.2 Comparison of radial velocity on the central axis of meridian plane

由图可以看出，所有曲线走向相同，近似正弦曲线，实验结果所得径向速度大于数值模拟结果。这是由于PIV法测得的是空间上的平均值，而实验中不可避免的存在振动和二次流的影响，因此在寻求与实验结果相似的数模结果时，所使用的雷诺数要大于临界雷诺数。且图中也可以看出，向外壁的流速要大于向内壁的流速，空间上向内流动的区域宽度要大于向外流动的区域。

2.3 雷诺数Re=124时流场模拟结果

2.3.1Re=124时流态分布

图3为雷诺数等于124时的子午面速度矢量图，上实线代表泰勒库特流装置内壁，下实线为外壁，内部曲线为速度等值线，等值线之间的差值为0.1倍内壁转速。

图3-a为实验结果、图3-b为数值模拟结果。由图可以看出液体总的流动走向与临界雷诺数时结果相似，近内壁处的高流速流体向外壁流动，近外壁处的较低流速的流体向内壁方向流动。且雷诺数为124时，图中泰勒涡明显增强，涡状流动更为显著。径向流速越大，等值线变化幅度越大，故两个相邻的泰勒涡之间等值线变化幅度最大。

2.3.2Re=124时径向流动特性

图4为雷诺数Re=124时，子午面中轴线上径向速度对比图。图中模拟结果和Wereley实验的结果拟合较好。

由图4所示，不论在走向上还是数值大小上，数值模拟结果与实验结果都较为相似。雷诺数Re=124时，其模拟结果仍满足由内壁向外壁流动径向流速最大，但该区域较小。由外壁向内壁流动区域流速较小，但该区域明显大于向外流动区域。这一结果符合质量守恒定律，向外壁的流速较大则可能是受到了离心力的影响。当雷诺数Re=124时，最大径向流速约是雷诺数Re=105时的三倍，即随着雷诺数的增大，子午面径向流速随之增大。

3 水平同轴转柱体间隙区域间流场数值模拟

3.1 水平同轴转柱体模型建立

为便于今后的制作考虑，本文数值模拟的水平同轴转柱体外圆柱直径设为184.1 mm，内圆柱直径设为165.1 mm，间隙宽度设为9.5 mm，轴向长度设为400 mm。对水平放置的同轴转柱体进行数值模拟，建立新的数学模型，在径向、轴向和周向上网格数量确定为32个、128个和256个，内壁面、外壁面和两个端面共四个面设置膨胀层，近壁面起始网格高度设置为0.1 mm，膨胀系数为1.1。

线型稳定性分析和能量理论常被用来计算临界泰勒数[19]，由此计算泰勒库特流装置的临界泰勒数，进而转化为临界雷诺数。经计算该装置由层流转变为层流泰勒涡时的临界雷诺数为124，流体完全发展为湍流的临界雷诺数为3 844，并以此来控制流态变化的模拟。

3.2 库特泊肃叶流模拟结果

当雷诺数低于临界雷诺数时，装置中的流场表现为层流状态，此时称为库特泊肃叶流。本文数值模拟雷诺数Re=65时同轴圆柱间的流场，模拟中将内圆柱的转速设置为0.082 7 rad/s，进行稳态层流的数值模拟，模拟结果子午面速度云图如图5。

由图可见，当内圆柱的转速较小，即所设雷诺数低于临界雷诺数时，整个流场简单而稳定，除了周向，其他方向上的流速基本为零，流体没有出现扰动，所产生的影响可忽略不计。此时同轴圆柱间隙间的流场呈现出典型的层流特征，满足库特泊肃叶流的要求。

3.3 层流泰勒涡模拟结果

将内圆柱转速提高到0.191 7 rad/s，此时的雷诺数为150，大于所计算出的临界雷诺数。间隙内的流动状态随之发生改变，流体受到扰动影响，由库特泊肃叶流发展为完全的层流泰勒涡。

3.3.1 层流泰勒涡流动特性

当雷诺数增加至150时，内外圆柱间隙内流场子午面速度云图结果如图6所示。

此时的流态发生了变化，整体出现了规则的扰动，速度等值线发生了规则的弯曲变形。由图可以看出，流场内流速仍有明显的分层，且较为均衡稳定。间隙内流场表现出层流泰勒涡的所有特征，图6直观地展示了层流泰勒涡在间隙内的流动形式。层流泰勒涡在运动时，流场内出现了规则分布的漩涡，但在此雷诺数下，涡状流速较小，周向流速仍为主要流动。

3.3.2 层流泰勒涡涡场

图7为雷诺数为150时子午面速度矢量图。上、下两条直线分别代表外壁和内壁，轴向上的曲线为速度等值线，每两条等值线的差值均为0.1倍的内圆柱转速。在轴向方向上可以看出有三个完整的涡状，流速等值线随着流动而发生变形。

如图所示，相邻的两个泰勒涡间速度等值线变形最大，即相邻的两个泰勒涡间径向流速最大，轴向流速最小。向外壁的最大径向速度要大于向内壁的径向速度，根据质量守恒定律，向内壁流动的区域大于向外壁流动的区域，即向内壁径向流动的两个泰勒涡之间的间距要大于向外壁径向流动的涡间的间距。

3.4 湍流泰勒涡模拟结果

进一步增大内圆柱的转速，当内圆柱转速增至6.390 7 rad/s时，雷诺数增加至5 000，此时间隙内的流动形态为湍流泰勒涡。

3.4.1 湍流泰勒涡流动特性

图8为流态变为湍流泰勒涡时子午面速度云图。由图可见，同轴圆柱间隙内流体出现了规则的扰动，速度更加均匀，不再具有明显的分层。

子午面速度云图直观地展示出间隙内湍流泰勒涡的流动形式。当同轴圆柱间隙内的主流变化为湍流后，轴向上仍然出现了规则分布的泰勒涡，与较低雷诺数下的子午面速度云图相比，涡的长度明显变长，间隙内泰勒涡的总数变少。

3.4.2 湍流泰勒涡涡场

图9为雷诺数为5 000时子午面速度矢量图，上、下两条直线分别代表外壁和内壁，轴向上的曲线为速度等值线，每两条等值线的差值均为0.1倍的内圆柱转速。图中包含三个完整涡的轴向部分，流速等值线随着流动变形明显。

图9 雷诺数Re=5 000时子午面速度矢量图Fig.9 Velocity vector of meridian plane at Re=5 000

如图所示，径向速度的最大值出现在相邻的两个泰勒涡间的交界区域，在一个泰勒涡的两侧的交界区域上径向速度方向相反。径向速度的最小值出现在通过涡心的径向线上，其值基本为零，而轴向速度则在此处达到了极值。与较低雷诺数下的子午面速度矢量图相比，二次流发展更加充分，周向速度更加均匀，在间隙内中心线附近区域周向速度大小几乎相同，在不同轴向位置径向线上，轴向速度的零点出现在不同的半径位置上，这说明了在湍流泰勒涡下，紊动更加剧烈。

3.4.3 湍动特性

湍动特性是湍流泰勒涡最基本的特性，本文采用单位质量内的湍动动能kt来描述湍动的大小。

(3)

图10为不同半径位置上的湍动动能沿轴向变化图，轴向上选取170～215 mm段，该段包含了一对泰勒涡，径向上选取的半径为84 mm、86 mm、87 mm、88 mm和90 mm位置处进行对比。

图10 不同半径位置上湍动动能沿轴向变化图Fig.10 Axial variation of turbulent kinetic energy at different radii

轴向上190～200 mm处，径向流动向外，内壁区域的湍动动能达到最大，轴向上170～175 mm和210～215 mm处，径向流动向内，外壁区域的湍动动能达到最大，且内壁附近的最大湍动动能要高于外壁附近的最大湍动动能。在径向流动向外的区域，半径越大的位置，湍动动能越小，而径向流动向内的区域与之相反，半径越大的位置处湍动动能越大。在泰勒涡中心部分，除了外壁附近区域，沿径向位置的湍动动能趋于一致，都达到最小值。

本节模拟水平放置的同轴转柱体间隙区域间流体流动，模拟较低雷诺数下的库特泊肃叶流与层流泰勒涡两种典型流场，以及高雷诺数下的湍流泰勒涡流场，库特泊肃叶流为周向层流，流场简单且较为稳定；随着雷诺数的增大，层流泰勒涡阶段，周向速度在径向上仍具有明显分层，同轴圆柱间隙内在轴向上出现规则分布的泰勒涡；在湍流泰勒涡阶段紊动剧烈，流速更为均匀，不再分层明显，泰勒涡中心部分，湍动动能基本一致。从总体上看，竖直与水平放置下的流场中，流态变化过程基本一致，模拟结果与叶立[20]等对泰勒库特流流态过渡变化的数值模拟结果相符。

4 紊动作用下的絮凝研究讨论

水流动力环境是细颗粒泥沙絮凝的关键影响因素之一，水流的紊动作用能够带动和加强泥沙颗粒间的碰撞从而促进絮凝的发生；另一方面，水流的强紊动作用会破坏絮团结构从而抑制絮凝的发生[21]。水流紊动可被看成各种不同尺度的涡旋运动与平均流速的叠加，涡旋运动产生的剪切力和惯性离心力是泥沙颗粒絮凝发生和破坏的主要作用力[22]。泰勒库特流装置具有剪切力小、对絮体破坏小、流动与停留时间易调等优点，因而常用于絮凝实验研究。

柴朝晖[4]利用同轴旋转双筒以层流转变为层流泰勒涡时的临界速度ω1、流体完全发展为湍流时的临界速度ω2为分界点，划分流态及水流强度。实验结果证明，当水流强度较小、流态以层流为主时，泥沙颗粒运动以絮凝为主；当水流强度较大、水流为湍流时，紊动作用强，阻碍泥沙的絮凝。韩晓婷[23]数值模拟泰勒库特流装置并做絮凝试验，结果表明泰勒涡涡旋结构清晰时絮凝效果最佳。内筒转速较小时，湍动动能较小，不足以使絮体颗粒发生有效碰撞，随着转速的增加，过大的湍动动能则会破坏絮体结构。

本文的数值模拟结果如表1所示，随着内圆柱转速的增加，同轴转柱体内流速、湍动动能逐级增加，流态随之变化，经历了由层流到湍流泰勒涡的过渡过程，本文模拟结果整体趋势与韩晓婷[23]数值模拟结果相似。竖直放置的同轴转柱体不可避免的存在泥沙浓度梯度，浓度梯度的存在必然会对絮凝产生一定的影响。本文提出将同轴转柱水平放置，利用周向流消除浓度梯度影响，更直观研究流态、水流强度的变化对絮凝的影响，因而更适用于絮凝实验，当然具体的结果还需要后续根据本文数值模拟结果制作物理模型并开展进一步论证研究。

表1 不同转速下数值模拟结果Tab.1 Numerical simulation results at different speeds

5 结论

(1)数值模拟结果与实验结果对比，数值模拟得到的结果与实验所得基本相同，两者差约为0.9%。考虑到误差存在的必然性，可以认为本文所建立的数学模型能较好地模拟出泰勒库特流装置内的流场特性，其结果较为精确。

(2)泰勒库特涡流场中，内圆柱附近流体流速较高，外圆柱附近流速较低，高速流体向外圆柱流动，低速流体向内圆柱流动，从而形成涡状。子午面径向流速最大值出现在由内圆柱向外圆柱流动的区域。液体由内向外流动的区域较小，由外向内流动的区域较大。

(3)对水平的同轴转柱间隙区域间流场进行模拟结果表明，随着内圆柱转速的提高，雷诺数逐渐提高，当雷诺数大于临界雷诺数时，同轴圆柱间的流动不再是稳定的层流，在轴向上会出现明显的均匀分布的泰勒涡，这一现象与竖直放置的泰勒库特流装置结果是相似的。

(4)随着雷诺数的增加，泰勒涡发展逐渐充分，涡状特征逐渐明显，水流流态也由层流变为涡流。雷诺数越大，水流涡动越强，径向流速越大，对流体的影响更为明显。但不同雷诺数下，径向流速与转速相比只占2%～7%，流场中转向流速仍起主导作用。

(5)径向流动向外的涡间区域，内圆柱附近湍动动能最大，沿径向向外半径越大的位置，湍动动能越小；与之相反，径向流动向内的涡间区域，外圆柱附近区域湍动动能最大，湍动动能由外向内逐渐减小；涡中心处，湍动动能基本一致。

(6)同轴转柱研究水流紊动对细颗粒泥沙絮凝的影响，流态为层流时，水流强度较小，泥沙絮凝；流态为湍流时，水流强度较大，泥沙絮凝受阻。水平放置的同轴转体，能有效避免竖直方向上的泥沙浓度梯度对絮凝的影响，可进一步开展论证研究，为今后水平同轴转柱体利用于絮凝实验提供支撑基础。