引言

静态混合器是一种依靠管道内周期排列的螺旋元件对流体不断的分割、变形、移位及混合的在线、高效混合设备，与搅拌混合器相比，具有结构简单、空间占有率小、停留时间短、耗能低及维修成本低等优点^[1-2]，因此广泛应用于精细化工、水处理、制药工业、生物化工等工业过程^[3]。

掌握静态混合器内流场的湍流特性规律对现有静态混合器的结构优化及新型静态混合元件的开发具有重要意义^[4-6]。Tajima等^[7]采用Kenics型静态混合器(KSM)对水和液态CO₂生成CO₂水合物的混合过程进行实验研究，结果表明湍流状态下随着KSM元件数量增加，液滴尺寸逐渐变小。Jaworski等^[8]采用大涡模拟(LES)和非定常雷诺平均NS方程(RANS)两种方法，对KSM内两相流的湍动能和耗散率等参数估计进行了有效性分析。Regner等^[9]利用计算流体动力学(CFD)软件分析了层流状态下不同长径比的Lightnin静态混合器(LSM)内的流场结构，研究发现带有曲率的混合元件交错处会产生径向二次流旋涡。本课题组前期采用大涡模拟分析KSM内二次流径向涡的形成诱因及演化过程，从多纵向涡的角度探讨了扭旋叶片对流体力学性能及传热性能的影响机理^[10]；并对LSM中高黏度流体的层流流动和混合特性进行数值模拟，发现LSM内高黏度流体的混合性能优于KSM^[11]；随后分析了LSM内湍流的传热特性和旋涡演化过程，结果表明半径尺度的强制涡受到LSM混合元件的影响逐渐耗散且被诱导形成两对迪恩涡，强化了LSM的传热性能^[6]。

本征正交分解(POD)方法被许多学者用于对湍流流场的测量结果进行分析，Rempfer等^[12]将POD应用于平板边界层内的相干结构和动力学行为的研究中，结果表明高阶的流动相干结构对转捩的产生有着十分重要的作用；余颖等^[13]利用POD方法提取内置蜗杆传热管流场中的拟序结构，发现多对大尺度纵向涡结构的产生会提高强化传热效率；张建伟等^[14]对双组分层对置撞击流的湍流速度场进行POD分解，发现低阶模态能够准确还原流场中的涡结构。

到目前为止，关于LSM的研究主要集中在湍流传热、层流混合等方面，而针对LSM的瞬态湍流特性及分布混合特性的研究尚未见文献报道。因此，本文利用POD方法提取流场的湍流结构，以期明晰LSM内流体的流动特性。采用平均涡量、拉伸率等参数表征LSM内流体的混合特性，为丰富LSM混合机理提供理论依据，并扩大其在工业上的应用范围。

1 数值模型 1.1 控制方程

大涡模拟是流体力学中一种重要的数值模拟方法，其基本思想是直接计算大尺度涡，对于小涡则通过亚格子模型进行模拟。区别于直接数值模拟(DNS)和RANS方法，LES在对大涡解析的同时还对小涡进行模化，通过空间滤波操作过滤掉小涡，滤波后大涡控制方程为^[15]

$ \begin{aligned} \frac{\partial\left(\boldsymbol{u}_{i}\right)}{\partial t}+\frac{\partial}{\partial x_{j}}\left(\boldsymbol{u}_{i} \boldsymbol{u}_{j}\right)=-\frac{1}{\rho} \frac{\partial p}{\partial x_{i}}+\nu \nabla^{2} \boldsymbol{u}_{i}-\frac{\partial \boldsymbol{\tau}_{i j}}{\partial x_{j}} \end{aligned} $

(1)

$ \frac{\partial \boldsymbol{u}_{i}}{\partial x_{i}}=0 $

(2)

式中，x为笛卡尔坐标系；u_i和u_j为流体速度；ρ、p和ν分别为流体的密度、压力和运动黏度。

亚格子应力τ_ij定义为

$\boldsymbol{\tau}_{i j}=-2 \mu_{t} \boldsymbol{S}_{i j}+\frac{1}{3} \delta_{i j} \boldsymbol{\tau}_{k k} $

(3)

其中，

$ \boldsymbol{S}_{i j}=\frac{1}{2}\left(\frac{\partial \boldsymbol{u}_{i}}{\partial x_{j}}+\frac{\partial \boldsymbol{u}_{j}}{\partial x_{i}}\right) $

(4)

式中，S_ij为应变速率张量，τ_kk为亚格子应力的各向同性部分。Smagorinsky-Lilly模型定义的亚格子湍流涡黏系数为^[16]

$ \mu_{t}=\rho L_{\mathrm{s}}^{2} \sqrt{2 \boldsymbol{S}_{i j} \boldsymbol{S}_{i j}} $

(5)

式中，L_s为亚网格尺度的混合长度。

1.2 本征正交分解

本征正交分解方法由Lumley^[17]首次提出，经Sirovich^[18]改进后发展为“快照POD”，逐渐被广大学者用于复杂流场的降阶及特征提取。假设在混合器横截面内有m个离散采样点，每次采样同时获取每个采样点的三维速度值u(x_m)。因此，获取N张快照脉动速度场数据，以矩阵U的形式表示为^[19]

$ \boldsymbol{U}=\left[\begin{array}{lllll} \boldsymbol{u}^{1} & \boldsymbol{u}^{2} & \boldsymbol{u}^{3} & \cdots & \boldsymbol{u}^{N} \end{array}\right]=\left[\begin{array}{cccc} u_{1}^{1} & u_{1}^{2} & \cdots & u_{1}^{N} \\ u_{2}^{1} & u_{2}^{2} & \cdots & u_{2}^{N} \\ \vdots & \vdots & \vdots & \vdots \\ u_{m}^{1} & u_{m}^{2} & \cdots & u_{m}^{N} \end{array}\right] $

(6)

式中，上标N表示快照数目，下标m代表三维速度的矢量位置。自相关矩阵为

$ \boldsymbol{C}=\boldsymbol{U}^{\mathrm{T}} \boldsymbol{U} $

(7)

速度相关矩阵的特征向量和特征值为

$ \boldsymbol{C A}=\lambda \boldsymbol{A} $

(8)

用式(8)计算得到的特征值及对应的特征向量去构建POD模态$ \phi^{i}$的基函数

$ \phi^{i}=\frac{\sum\limits_{n=1}^{N} \boldsymbol{A}_{n}^{i} \boldsymbol{u}^{n}}{\left\|\sum\limits_{n=1}^{N} \boldsymbol{A}_{n}^{i} \boldsymbol{u}^{n}\right\|} $

(9)

式中，i=1, …, N，A_nⁱ为特征值λ_n对应的特征向量。

POD系数

$ \boldsymbol{a}^{n}=\boldsymbol{\psi}^{\mathrm{T}} \boldsymbol{u}^{n} $

(10)

线性重构任意时刻的速度场

$\boldsymbol{u}^{n}=\boldsymbol{u}(x)+\sum\limits_{i=1}^{N} a_{i}^{n} \phi^{i} $

(11)

1.3 物理模型

假设圆管直径D=40 mm，总长度为560 mm，对含有元件长径比A_r分别为1.0、1.5、2.0和3.0的LSM流场进行比较分析。LSM内混合元件由两个半椭圆形板和两个三角形板组成，相邻两元件首尾相连且相互交错90°布置安装，几何模型如图 1所示。对于长径比A_r=1.5的LSM，两个交叉半椭圆板分别与轴截面成34°和146°，其他模型的具体几何参数见表 1。

1.4 模拟策略

本文采用ANSYS FLUENT V16.1软件中的LES湍流模型，对图 1中4种不同长径比的LSM进行非稳态数值模拟。以水(密度ρ=998.2 kg/m³，黏度μ=0.001 003 Pa·s)作为工质。采用课题组前期对LSM湍流流动传热预测中的网格无关性验证结果^[6]，利用多面体网格对计算域进行离散化处理，采用10层边界层和1.2的网格增长率对壁面进行细化，总网格数量为1 597 947。

采用三维压力基求解器，并使用SIMPLEC算法对压力和速度耦合；梯度项求解采用Least Squares Cell Based离散法；压力离散采用标准形式，动量离散采用Bounded Central Differencing；采用用户自定义函数(user define function，UDF)设置速度入口边界^[6]，对应的雷诺数分别为1 990、3 981、5 971和7 962。出口边界条件设置为outflow。连续性方程收敛残差小于10^-4，其余变量的收敛残差为10^-6。

1.5 模型有效性

大涡模拟在降低空间和时间分辨率的前提下可以对具有较高雷诺数和较复杂的湍流流动进行模拟，计算精度高。课题组前期采用LES方法对KSM内的流场进行研究，数值计算的速度分布和实验一致性较好^[10]。LSM和KSM混合器都属于螺旋流静态混合器^[9]，具有相似的工作原理。因此，本文选用大涡模拟对Lightnin型静态混合器进行研究。

2 POD分析与结果讨论 2.1 采样数量选择

为探究采样数量对POD分析精确度的影响，分析了不同采样数量下的模态能量分布，如图 2所示。从图中可以看出，随着模态阶数s的增加，POD模态的特征值λ迅速衰减，当模态阶数超过10之后，高阶模态所携带的能量较低。其中，前四阶模态携带的能量最多；在采样数量超过2 000时，前四阶的模态能量占整体的29%以上，且随着样本数量的增加低阶模态所包含的能量不再发生变化。因此，综合考虑计算精度及成本，选取采样数量为2 000。

2.2 POD模态收敛性验证

对LSM数值模拟所得的样本进行计算，评估不同模态下的二范数ε_L2(L₂-norm)^[20]，具体公式如下。

$ \varepsilon_{L 2}=\left\|\delta_{i, 1}-\delta_{i, 2}\right\|_{2} $

(12)

式中，δ_{i, 1}、δ_{i, 2}分别为两次不同采样频率样本中POD分解之后的第i阶模态。

样本数量与模态收敛性之间的关系如图 3所示。由图可知，模态之间的差异随着样本数量的增加逐渐减小。当样本数量n小于500时，差异值变化明显，在样本数量很小的情况下，L₂-norm可达到1.3左右。当n>1 000时，L₂-norm小于0.3，并且逐渐趋于0，表明经POD得到的模态差异性稳定，即POD模态是收敛的。

2.3 POD模态能量分析

图 4为不同轴截面LSM内湍流流体前四阶模态的能量分布。从图 4可以看出，不同阶的模态能量在z/l=0~2范围内呈减小趋势，表明流场复杂度不断增加。在流场复杂度稳定后，1阶模态能量所占的湍动能分数在8%~10%，2阶模态能量所占湍动能分数在6%~8%，3阶模态能量所占湍动能分数在5%~7%，4阶模态能量所占湍动能分数在5%~6%。前四阶模态能量占总湍动能的24%~31%，包含了主要的大尺度旋涡结构信息。

图 5为前四阶POD模态的径向速度矢量图。由于1阶模态包含着最多的流场信息，从图 5(a)可以清晰地看出在元件壁面存在两对大尺度旋涡结构且对称分布，周围无小尺度旋涡结构。随着模态能量的降低，在第2、3阶模态中小尺度旋涡数量逐渐增加，所提取的流场信息逐渐减少，这主要是由流场的复杂程度增加所致。在第4阶模态中旋涡结构开始散乱，小尺度的旋涡结构尺寸逐渐增大。

对LSM内的原始瞬态流场(图 6(a))进行40阶本征正交分解，利用前四阶模态重构流场流线图，如图 6(b)所示。对比发现，原始流场在管壁边缘处存在小尺度涡对，干扰流场大尺度相干结构的提取与分析。POD方法在时均速度场的基础上对任意时刻的原始流场进行线性重构，能够过滤小尺度的旋涡结构，同时保留半圆形流道内原始流场的特征，使核心流区域的径向二次流大尺度相干结构的空间分布更加明晰。

2.4 回流系数

二次流能够强化静态混合器的混合性能。为了更好地了解LSM内的二次流分布，通过计算整个采样过程流场中选定位置处流速方向与主流方向相反的时间与总采样时间的比值得到回流系数γ_t。图 7为在z/l=4.1~6.0处LSM内的回流系数云图。由图中的高亮区域可以看出，回流主要发生在管壁与元件连接处附近以及元件中心截面处。

2.5 涡量分布

图 8描述了不同雷诺数下平均涡量〈ω〉随轴向位置的变化。从图中可知，受到元件旋向影响，不同轴截面内平均涡量分布具有正弦状的周期性变化规律。在混合器入口段无元件分布，没有旋涡产生，〈ω〉为0；在z/l=4.0~5.0范围内，〈ω〉逐渐增大，表明在一个混合元件周期内旋涡从正向涡逐渐向反向涡转变，峰值出现在左旋和右旋元件连接处；在z/l=5.0~6.0范围内，〈ω〉逐渐减小，表明旋涡从反向涡逐渐向正向涡转变。随着雷诺数的增加，涡量峰值出现不同程度的变化，当Re=1 990时，正向涡值的峰值平均值为6.37，而Re=3 981、5 971和7 962对应的正向涡值峰值平均值分别为15.83、25.33和36.39，涡量分别增大了148.5%、297.6%和322.7%，表明增大Re可提高LSM内流体的湍流程度。

2.6 混合特性

前期的研究成果表明，静态混合器内的流动具有混沌混合特性^[21]。采用拉伸率Λ_n这一指标对不同长径比和不同雷诺数下的LSM分布混合性能进行定量分析，拉伸率数值越大，代表强化分布混合能力越强。根据Ottino混合理论^[22]，每个粒子所经历的拉伸量可以通过小矢量的变形求出。

图 9为不同长径比和雷诺数下LSM内粒子拉伸率的几何平均值。由图 9可以发现，LSM内粒子的拉伸率均随轴向长度的增加而线性增长；在同一轴向位置，A_r=1.0对应的拉伸率最大，A_r=3.0对应的拉伸率最小，说明长径比越小，拉伸效果越明显。这是由于相同轴向距离内小长径比的混合叶片扭曲程度更大，径向二次流更明显，并且在相同长度的管路内，可以布置的小长径比的元件个数多于大长径比的元件个数，因此，流体微元在经过布置有小长径比元件的管路中时被拉伸的概率更大。从图 9(b)可以看出，随着雷诺数的增加，LSM的最大拉伸率和横截面的平均拉伸率均增大，Re=3 981、5 971和7 962下的拉伸率与Re=1 990时的相比，分别提高了11.5%、15.1%和18.9%。

通过对LSM内的拉伸率数值进行统计学分析，可以定量评价流体的微观混合程度^{[21, 23]}。图 10为LSM在不同长径比和不同雷诺数下出口截面的拉伸率概率密度分布。从图中可以看出，对流体有切割分流作用的LSM，其拉伸率概率密度分布曲线形状类似高斯分布曲线，说明拉伸率数值分布较为集中；径向二次流随着长径比的变小而增强，致使概率密度曲线向着高拉伸率方向移动，这是由于流体经过较大的拉伸变形，产生较高的拉伸率。A_r=1.0时曲线在尾部的震荡相较于其他长径比的要明显，这说明其内部的高拉伸率的粒子数量较多。随着长径比的增加，各曲线峰值逐渐向左移动，说明小长径比的混合元件对流体的混合能力较强。随着雷诺数的增加，拉伸率分布范围没有改变，但各个截面处平均拉伸率的概率密度增大，总体上提高了混合器的混合能力。

采用压降Δp和混合性能强化系数η来分析不同长径比的LSM的流动混合强化行为^[11]。不同长径比LSM内的压降及混合性能强化系数的变化如图 11所示。可以看出压降Δp和混合性能强化系数η均随长径比的增大逐渐减小。在Re=3 981条件下，A_r=1.0的LSM内的压降Δp=712.9 Pa，混合性能强化系数η=6.8。随着长径比增加，A_r=1.5、2.0、3.0的LSM的压降分别是A_r=1.0时的28.5%、14.2%和0.06%，而混合性能强化系数分别是A_r=1.0时的93.5%、85.2%和66.9%。综上所述，在A_r=1.5时LSM的压降明显减小，同时具有较好的混合性能。

3 结论

(1) 采用POD方法对LSM内的瞬态湍流流场进行分析，前四阶模态能量占总湍动能的24%~31%，蕴含主要的大尺度旋涡结构信息，可消除小尺度旋涡结构对流场的干扰。

(2) 对于同一长径比的LSM，涡量随着元件的轴向位置变化呈正弦状周期性变化，随着雷诺数的增加，不同轴向截面的平均涡量均增大。

(3) 对拉伸率的分析结果表明，LSM内部的拉伸变形随着雷诺数的增加而增加；流体微元在经过布置有小长径比元件的管路时被拉伸的概率更大。

(4) A_r=1.0时LSM的混合强化性能最好，然而与A_r =1.5、2.0、3.0时相比压降较高。综合压降和混合性能强化系数，A_r=1.5的LSM具有较好的混合性能，同时压降降低显著。