Analytical Study for Uncoordinated Deformation of Existing Pipeline and Soil Induced by Tunnel Undercrossing

Taylor-Couette涡流^[1]在不稳定现象及流动转捩的研究中扮演着重要的角色. 自Taylor以后，研究者从实验测量^[2-4]和数值模拟^[5-6]等不同角度对Taylor-Couette流进行了研究，且将其应用于黏度测量、絮凝、混合搅拌、旋转式反应与分离等多个领域. 实验测量面临着诸多技术限制. 一方面，实验室制作出的实验模型会受到制作技术、结构运行等方面的干扰，降低流场反应精度；另一方面，在测量过程中，环境噪声干扰也会放大测量误差，导致实验产生偏差. 随着科学技术的不断发展，用数值模拟研究流场特性可以弥补上述不足，也是十分必要.

Taylor-Couette反应器由2个同心转筒构成，若保持外筒固定，随着内筒转速增加，两筒环隙间流体会依次出现层流涡、波状涡、调制波状涡和湍流涡等不同的流态^[7]. 在前期的研究中，发现波状涡流场可以为絮凝反应提供最佳的水力条件^[3]，而且波状涡流态在其他领域也一直备受关注^[8]，有必要进行深入探索. 为更细致地了解Taylor-Couette波状涡流场，本文采用大涡模拟方法（large eddy simulation, LES）^[9-11]对Taylor-Couette反应器环隙内的波状涡流场进行瞬态数值模拟，从二维和三维的角度研究环隙流场的波动变化特征.

1.   数值模拟及分析方法

1.1   大涡模拟基本控制方程

大涡模拟方法综合了直接数值模拟和雷诺平均数值模拟的优点，将流体运动分为可解大尺度运动和不可解亚格子尺度运动，对大尺度可解部分的紊流通过求解方程进行计算，而用亚格子尺度模型来弥补亚格子尺度运动对大尺度运动的影响. 对于不可压缩流动，通过滤波函数对不可压缩Navier-Stokes方程进行过滤，得到大涡模拟控制方程^[12]为

式中：ρ为流体密度；t为时间；i、j为张量角标（取值为1、2、3）；${\bar {\boldsymbol{u}}_i}$、${\bar {\boldsymbol{u}}_j}$为流体分别沿坐标轴${x_i}$、${x_j}$方向经过过滤后得到的速度矢量；μ为流体动力黏性系数；$\bar p$为经过过滤后得到的大尺度压强；${{\boldsymbol{\tau}} _{ij}}$为亚格子尺度应力张量^[13]，是不封闭项，需对其进行建模封闭，采用Smagorinsky-Lilly亚格子模型，如式（4）所示.

式中：${\mu _{\mathrm{t}}}$为亚格子湍流黏度；C_s为Smagorinsky常数；$\varDelta$为滤波宽度；$\left| {{{\overline {\boldsymbol{S}}}}} \right| = \sqrt {2{{{{\overline {\boldsymbol{S}}}}}_{ij}}{{{{\overline {\boldsymbol{S}}}}}_{ij}}}$，为大尺度的应变值；${\overline {\boldsymbol{S}} _{ij}}$为经过滤波后的速度变形张量；k为亚格子湍动能；δ_ij为克罗内克尔算子.

1.2   几何模拟及数值方法

旋转雷诺数Re的改变直接影响Taylor-Couette流场形态，在环隙间液体的运动黏度一定的情况下，Re只与内筒转速n有关，如式（7）所示.

式中：r₁为内筒半径，d为环隙宽度，$\nu$为液体的运动黏度（20 ℃时水的运动黏度为1.006 × 10⁻⁶ m²/s）.

Taylor-Couette反应器几何模型及网格划分细节如图1所示. 内筒旋转，外筒固定，模型几何尺寸如下：内筒半径r₁=37.5 mm，外筒半径r₀=50.0 mm，环隙宽度d=r₀−r₁=12.5 mm，内外筒半径比η=r₁/r₀=0.75，筒高L=440 mm. 结合研究目标与对象的属性，设置第1层网格高度为0.02 mm，分别划分10层增长率为1.2的边界层网格. 在满足网格无关性的要求下，综合考虑网格质量、计算精度和计算速度，确定网格总数为209万个.

图  1  Taylor-Couette反应器几何模型及网格划分示意

Figure  1.  Geometry model and meshing for Taylor-Couette reactor

下载: 全尺寸图片 幻灯片

用Fluent 2020作为计算主体，模拟方法选择LES，亚格子尺度模型选择Smagorinsky-Lilly，C_s常数经数值模拟方法验证后最终确定为0.1；边界条件：上下底面为自由表面，内外筒壁面分别为旋转和固定边界；压力速度耦合选择SIMPLE格式，空间离散化梯度项选择基于单元体的最小二乘法，压力项为二阶格式，动量项为有界中心差分格式，瞬态离散方案选择二阶迎风格式.

控制内筒转速，对转速5～45 r/min（旋转雷诺数Re为244～2196）的Taylor-Couette流场进行瞬态数值模拟，即可以得到不同形态的涡流场. 为在有限的篇幅内阐明问题，在具有明显波状涡特征的流场中选取转速分别为10、20、30、40 r/min共4个工况下的模拟结果进行分析.

1.3   大涡模拟与实验测量验证

参考对照文献[3]中利用粒子图像速度场仪（particle image velocimetry, PIV）的测量方法，引用PIV测量流场数据与数值模拟结果对比. 图2为各转速下大涡模拟和PIV测量的环隙子午面速度矢量图对比. 由图可以看出：大涡模拟得到的矢量图与PIV测量结果在涡旋形态上较为相似，一定程度上可以反映大涡模拟结果的可靠性；此外，大涡模拟结果相对规整有序，而PIV测量矢量图则因为外界因素的干扰，显得略微杂乱，某些局部区域由于噪声甚至出现矢量紊乱现象. 综上所述，大涡模拟可以较为准确地反映环隙间流场的流动特征.

图  2  大涡模拟与PIV测量环隙子午面速度矢量图对比

Figure  2.  Comparison between LES and PIV measurement of velocity vector on meridian plane

下载: 全尺寸图片 幻灯片

1.4   涡流场涡对界定及特征值表征

涡对的界定便于追踪研究目标，并从多维角度对Taylor涡展开研究. 将环隙子午面上两外向流（由内筒流向外筒）之间形成的2个旋转方向相反的涡旋界定为一对涡对，如图3所示.

图  3  涡对及特征线位置示意

Figure  3.  Vortex pairs and position of characteristic lines

下载: 全尺寸图片 幻灯片

为详细表征Taylor-Couette波状涡流场的速度波动变化特征，选取矢量图中任意一涡对，在涡对内两涡交界位置处作特征线1，涡对间两涡交界位置处作特征线2. 在每条特征线上等距离选取10个特征点，从左至右（外壁到内壁）依次用字母A～J表示（特征点颜色与后续对应图的图例颜色一致）. 如图3右侧局部放大图所示，特征线1代表涡对内涡间的流动特征（内向流：由外筒流向内筒）；而特征线2代表涡对间的流动特征（外向流：由内筒流向外筒）. 图中只示意了各特征线与涡对位置的关系. 特征线在环隙的高度（轴向距离）是动态的，即不同工况下，特征线的坐标数值会有所不同，具体由两涡交界的位置决定. 随着转速或时间的变化，涡的大小和位置均会改变，导致特征线在坐标轴上的位置相应发生变化.

1.5   涡识别方法

依据刘超群^[14]的Ω涡识别方法，将总涡量分解为旋转部分涡量和非旋转部分涡量，引入参数Ω表示旋转部分涡量占总涡量的比例，不仅避免了人为调整阈值的问题，还能同时捕捉强涡和弱涡. Ω的表达式为

式中：a、b分别为对称张量、反对称张量，是基于速度梯度$\partial {{\boldsymbol{u}}_i}/\partial {x_j}$的特征值演化得来的，对称张量反映流体的变形效应，反对称张量反映流体的旋转效应；${\Vert {\text{•}}\Vert }_{\text{F}}$为矩阵的Frobenius范数；ε为极小正数，以避免当分母为极小数时出现极大误差^[15].

显然，Ω∈[0,1.0]，表示涡量的浓度. 当Ω=1.0时，流体做刚体旋转；当Ω>0.5时，表示反对称张量b相较于对称张量a占优.

2.   模拟结果

2.1   二维环隙子午面速度矢量场波动变化特征

图4为10、20、30、40 r/min工况下环隙子午面速度矢量图的周期变化过程展示. 不同工况下波状涡流场周期不同，随着转速增大，周期缩短，各工况周期时间用T_n进行区分. 为便于观察与比较，选取环隙子午面中任一涡对在一个周期内不同时刻（如图4（a）中的0、T₁₀/4、T₁₀/2、3T₁₀/4、T₁₀，余图同）的速度矢量场，并将涡对内两涡旋分别命名为 ①、② 号涡旋，如图4所示. 图中：虚线箭头为 ①、② 号涡旋在一个周期中的变化过程指示线；蓝色实线箭头为不同时刻环隙间主流液体流动方向引导线. 可以发现，速度矢量图上涡旋大小、涡心位置和主流液体流动方向等物理量的变化情况存在周期性变化规律，在经历一个周期的复杂变化后，各工况下各物理量特征在周期始末时刻基本一致.

图  4  各转速工况下速度矢量图周期变化过程

Figure  4.  Periodic variations in velocity vector diagram at different rotational velocities

下载: 全尺寸图片 幻灯片

2.1.1   涡旋大小

如果将 ① 号涡旋面积最大时刻（对应 ② 号涡旋面积最小）视作周期开始，则周期起始阶段，图4中 ① 号涡旋的面积最大，② 号涡旋的面积最小；随着时间行进，① 号涡旋的面积由最大减至最小再增至最大，② 号涡旋的面积同步由最小增至最大再减至最小；周期结束时刻，2个涡旋的面积大小均又变回与周期起始时刻一致.

2.1.2   涡心位置

随着周期行进，图4中两相邻涡旋的涡心位置在竖直方向和水平方向上发生振动和漂移（竖直方向为振动，水平方向为漂移），并在周期结束时刻回到原点. 且相比于10、20 r/min的变化情况，30、40 r/min时的涡心位置变化幅度较为轻微一些. 到45 r/min以后，在一个周期内，两相邻涡旋的涡心位置在肉眼视觉上几乎观察不到其随时间发生的改变.

2.1.3   主流液体流动方向

当转速为10、20 r/min时，图4中涡间流体的流动方向发生了明显的周期性变化（蓝线所示），以转速20 r/min为例：初始时刻的流动方向为绕 ① 号涡旋逆时针向下流动（或是绕 ② 号涡旋顺时针向下流动）；随着周期行进，流动方向变为绕 ① 号涡旋逆时针向上流动（或绕 ② 号涡旋顺时针向上流动）；在周期结束时刻，涡间流体的流动方向又变回到与初始时刻流动方向一致.

通过上述分析可知，Taylor-Couette波状涡的速度矢量场具有周期特征，且不同工况下周期波动特征不相同. 综合分析各物理量的周期性变化规律发现，涡间流体流动方向的周期性变化导致其他物理量发生变化. 以转速为20 r/min时 ① 号涡旋（② 号涡旋变化趋势相反）为例：当涡间流体流动方向趋势向上时，涡心位置向上移动，涡旋变小；当涡间流体流动方向趋势向下时，涡心位置向下移动，涡旋变大. 当转速为30、40 r/min时，涡间流体流动方向随周期变化的波动幅度逐渐减弱，涡旋大小与涡心位置的波动变化也相应减弱了许多.

2.2   环隙子午面特征线上各向速度波动变化特征

为阐明Taylor-Couette波状涡流场的波动特征，以目标工况特征线上的各向速度为对象，选取一个周期内近外壁面的点A、B，环隙中部的点E、F和近内壁面的点I、J的各向速度变化情况进行对比分析，各工况对应周期均分为20个时刻值.

2.2.1   轴向速度波动变化

各工况特征线上的一个周期内，轴向速度随时间的变化情况如图5所示. 由图可知，各工况特征线上的轴向速度方向在一个周期内均会经历两次相反的变化，并在周期结束时刻变回到与初始时刻方向一致，体现在数值属性上为正-负-正或者负-正-负. 说明涡间流体的轴向速度方向在不断地发生改变，并具有周期波动特征. 而随着转速增大，速度数值的最大值均明显增大，特征线2的数值明显大于特征线1，且波动也更为剧烈，说明主流液体外向流动增强.

图  5  各转速工况特征线上一个周期内轴向速度随时间的变化

Figure  5.  Variation of axial velocity with time along characteristic line at different rotational velocities within a cycle

下载: 全尺寸图片 幻灯片

2.2.2   径向速度波动变化

各工况特征线上的一个周期内，径向速度随时间的变化情况如图6所示. 由图可知，各工况特征线上的径向速度方向在各自周期内均保持不变：特征线1上的径向速度方向均保持正（内向流），速度值也在较小幅度内波动，且随转速增大而增大；特征线2上的径向速度方向均保持负（外向流），方向与特征线1相反，相应点的速度数值比特征线1上的更大，且随转速增大而增大，其波动也更为剧烈. 说明环隙间波状涡流体外向流的趋势大于内向流.

图  6  各转速工况特征线上一个周期内径向速度随时间的变化

Figure  6.  Variation of radial velocity with time along characteristic line at different rotational velocities within a cycle

下载: 全尺寸图片 幻灯片

2.2.3   切向速度波动变化

各工况特征线上的一个周期内，切向速度随时间的变化情况如图7所示. 由图可知：各工况特征线上的切向速度方向在各自周期内均保持不变，其方向均为负（环隙间流体绕内筒逆时针旋转），速度值也在较小幅度内波动，随着转速增大其绝对值也顺序增大，且变化过程中各点速度数值互不超越；各点的切向速度绝对值由内壁到外壁（点J至点A）逐渐减小，外壁处点A在周期内各时刻的数值均远小于其他点，不过特征线2的数值明显大于特征线1. 说明环隙间各点切向速度大小主要由距内筒的距离决定，距离越远，切向速度值越小；另外还能说明环隙间流体切向流动的能量来源于旋转内筒.

图  7  各转速工况特征线上一个周期内切向速度随时间的变化

Figure  7.  Variation of tangential velocity with time along characteristic line at different rotational velocities within a cycle

下载: 全尺寸图片 幻灯片

2.2.4   各向最大速度值波动情况对比分析

为更直观地了解环隙间流体在涡间各个方向上波动传递的强弱情况，提取各工况特征线上一个周期内各向最大速度分别进行对比分析，见图8所示.

图  8  各工况特征线上各向最大速度对比

Figure  8.  Comparison of maximum velocities along characteristic line under various working conditions

下载: 全尺寸图片 幻灯片

由图8（a）可以发现：随着转速增大，轴向最大速度的绝对值均增大，说明流场中流体轴向运动趋势增强；并且特征线2上的绝对值大于特征线1，说明流体轴向运动在外向流的涡对间较强，在内向流的涡对内相对较弱. 由图8（b）可以发现：特征线2上的径向最大速度值大于特征线1，说明涡间流体外向流趋势大于内向流，且以外向流为主. 可理解为在外向流体触碰到外壁时，流体改变方向冲向内壁，之后当流体触碰到内壁时再次改变方向，这种碰撞促使流体发生旋转，并在两外向流之间产生2个旋转方向相反的涡旋，形成一个稳定的涡对. 由图8（c）可以发现：随着转速增大，特征线上的切向最大速度也增大，且特征线2上的切向最大速度绝对值均大于特征线1，同样说明流体在外向流的涡间流动更为剧烈.

综上所述，切向速度数值远大于轴向和径向速度，所以一定程度上说，流场中切向上的拖拽力相对较大，带动环隙流体绕内筒做旋转运动，同时在轴向与径向速度的耦合影响下，在环隙子午面上形成涡旋，且流体会在环隙空间沿轴向进行上下流动与传递. 另外，各特征线上各向最大速度均与内筒转速正相关，涡间各向速度值均随着内筒转速的升高而变大，波动幅度也增大.

2.3   三维环隙涡流场波动变化特征

上述内容分析描述了波状涡流场在环隙子午面上随时间和空间的周期性变化过程. 由于环隙空间中任何径向位置都存在子午面，从数学角度上说子午面数量是无穷的，所以从全流场空间视角可以更全面地分析波状涡的波动特征. 采用Ω涡识别方法对环隙间的三维速度矢量场进行涡识别. 令Ω=0.52，得到环隙空间的三维涡旋等值面图. 为便于分析各工况下三维涡旋等值面图的动态变化过程，选取一个周期内的变化情况进行展示，从三维角度对流场波动特征展开全面分析.

图9为各工况一个周期内环隙间全流场三维涡旋等值面图的波动变化过程. 各子图左侧系列图为环隙全流场一个周期内不同时刻的三维涡旋等值面图；右侧系列为对应工况下某一时刻的涡旋等值面局部放大图. 连续观察左侧系列图中不同时刻的变化情况，可以看出，环隙空间的涡旋等值面形成了以涡旋为单元的串联涡旋通道，涡旋通道再两两组合配对，形成互偶涡对（涡旋区域形态特征互偶），再以互偶涡对为单元形成立体螺旋偶合涡结构. 涡对内涡间的螺旋偶合特征明显，一个上凹一个下凸，一个收缩一个伸展. 根据涡旋等值面图中各物理量的变化特征可以得到各工况（10、20、30、40 r/min）的变化周期分别为12.94、6.80、1.93、1.49 s.

图  9  各工况三维涡旋等值面周期内变化

Figure  9.  Variation of three-dimensional vortex isosurface within a cycle under various working conditions

下载: 全尺寸图片 幻灯片

另外，不同转速下涡对的波动偶合幅度也不相同. 10、20 r/min的涡旋表面光滑，流场波动偶合特征明显；而30、40 r/min涡旋表面光滑度下降，难以捕捉到明显的波动偶合特征，且外向流涡对间的涡旋表面开始长出“牙点”状突起（图9中黄色虚线所示区域），甚至出现接触连接到一起的趋势. 对不同内筒转速下波状涡环隙全流场中的涡对数量进行统计发现：15、18、20、25 r/min时为15对，12、28 r/min时为16对，10、30、35、40 r/min时为17对. 可以看出，波状涡的涡对数量随内筒转速增大呈现先减少再增加的趋势. 说明Taylor-Couette波状涡流场中的涡对数量与涡形态密切相关，不同内筒转速下的涡形态不同，其涡对数量也不一样. 所以，从有限空间的角度上来说，环隙内部涡对多，应缘于其涡旋结构与涡间距小，才有生成更多涡对的可能；而如果涡旋结构与涡间距大，有限空间内不足以容纳更多涡对，环隙内部涡对数量相对要少.

综上，在波状涡范围内，当内筒转速较低或较高时，涡对结构较小，环隙间涡对数量多；当转速为波状涡中间区域时，涡对结构大，环隙间涡对数量少.

3.   分析与讨论

为更清晰地分析说明流体在三维环隙间的波动情况，以20 r/min（旋转雷诺数Re为976）时环隙子午面上的流动情况为例，将半周期内流体在涡间三维流动情况示于图10. 局部放大图中，红、蓝色曲线为数值模拟生成的不同时刻的涡间主流液体的流动流线. 由图可以看出，任意时刻t到半周期时刻T_n/2，流线的方向发生了改变，且在经历一个完整周期T_n后又恢复初始方向. 说明Taylor-Couette波状涡流场中存在主流液体的周期性流动，流动方向随时间发生周期性波动与改变，不同时刻下主流液体的流动会引起环隙间流体微团随涡旋向上或向下移动，微团在随子午面上的涡旋一起上下移动的同时也被主流液体裹挟进入到不同的涡旋中. 将这种运动趋势扩展到环隙三维空间中：涡旋偶合连接到一起在环隙空间中形成涡旋串联通道，通道中的涡旋不断上下移动促成涡旋通道做螺旋偶合旋转运动. 流体微团随螺旋偶合涡串联通道可以被裹挟进入到不同空间位置的任一涡旋中.

图  10  环隙子午面半周期内流线波动三维示意

Figure  10.  Three-dimensional diagram of streamline fluctuations in half-period on meridian plane

下载: 全尺寸图片 幻灯片

结合2.2节对涡间流体各向速度周期内波动变化特征分析可知，波状涡流场中的流体运动主要发生在外向流的涡对间，而内向流的涡间流体运动相对较弱，这些特征均与三维波动变化一致. 如图10所示（右侧图中以灰色箭头方向代表流动方向，箭头大小代表流动强弱），涡对a、b间的流体运动强于涡对内的两涡间. 另外，由图10还可以发现，涡对a、b间的间距小于涡对内两涡的间距；且涡间距越小，涡间速度（液体流动）变化越剧烈，在同等时间内流体微团被裹挟经历的涡数量越多.

4.   结　论

1） Taylor-Couette波状涡流场二维环隙子午面速度矢量场存在周期性波动变化特征，轴向速度方向在周期内至少发生2次变化；径向与切向速度方向在各自周期内均保持不变，但涡两侧的径向速度方向相反，流场内切向速度方向相同且由内壁到外壁逐渐减小. 随着旋转雷诺数增大，各向速度值与波动幅度也增大. 所以，轴向与径向速度的不断变化在子午面上促成涡旋形成，叠加上切向速度后促使环隙间流体微团绕内筒做螺旋偶合涡旋转运动.

2） Taylor-Couette波状涡流场二维环隙子午面特征线2上各向最大速度比特征线1的大，说明涡间流体以外向流为主，主流液体传递主要发生在外向流的涡对间，动量输运也主要在外向流区域.

3） Taylor-Couette波状涡流场三维涡旋偶合波动现象明显，具有周期特征，不同转速下涡对的波动幅度不相同，在波状涡的中间区域，涡对的波动幅度最大. 随着旋转雷诺数增大，波动周期缩短；各工况（10、20、30、40 r/min）的周期分别为12.94、6.80、1.93、1.49 s.