徐百平，刘彪，刘尧，谭寿再，杜遥雪，刘春太

（1 五邑大学智能制造学部，广东江门529020； 2 广东轻工职业技术学院广东省高分子材料先进加工工程技术研究中心，广东广州510300； 3 郑州大学橡塑模具国家工程研究中心，河南郑州450002）

引 言

同向双螺杆挤出机具有混合性能好、自清洁效果好、脱气效率高等优点，在聚合物加工、反应挤出、食品加工以及制药等领域获得了广泛应用。这类机器内，一对形状相同的内部转子，安置在水平腔内同向等速旋转，其截面由几段不同半径的圆弧组成。为了实现自洁功能，避免运动干涉，两个转子之间有特定的错列角。一对三角形转子类似于一对三头双螺杆元件和捏合块，在一定程度上与挤出机的捏合盘具有相同的混合原理。这样的几何结构，对于螺纹元件的横截面内混合也有一定的参考意义。然而，当需要更大的工作空间和热交换面积时，最简单的方法是缩小转子尺寸以提供更多的空间，但是通常会损失部分自清洁特性。例如，当选择一对三角形转子时，尽管失去了自清洗功能，但可以提供更多的空间和传热表面。而且当转子尺寸缩小后，可以改变错列角来完成不同的混合过程。但这样的尺寸放大缩小以及错列角的选择对混合室内的混合效率影响效果还需要进一步研究。

1997 年，Avalosse 等[1]通过实验观察和数值模拟研究了三角形转子对作用下的混合问题。随后Bertrand 等[2]提出了另一种虚拟内嵌边界方法，使用拉格朗日乘子法来施加移动边界条件，进一步进行局部自适应网格细化来精确尖角处边界条件，从而提高了混合演化细节预测。当前，计算流体动力学（CFD）取得了巨大的进步，为评估和优化反应器设计提供了一种很有效的方法。有限体积法（FVM）、有限元法（FEM）、光滑质点流体动力学法（SPH）等数值方法现已用于模拟混合装置中的复杂流动[3-20]。其中，基于有限元的网格叠加技术（MST）取得了巨大的进步，该方法针对流动域和内部运动部件分别进行网格划分，然后根据运动部件的运动规律，在每个时间节点更新运动部件的运动位置。Kokini团队[13-16]将有限元数值模拟方法和叠加网格技术相结合，开展了食品加工方面的混合机理研究，取得了丰硕的成果。

针对三角形转子，虽然前人的研究在描述混合细节和提高数值精度方面取得了积极进展，但对混合的机理理解仍然不够清晰。没有采用动力学研究工具，如统计分析、Lyapunov 指数及Poincaré 截面等手段，也没有探索转子大小和几何形状、错列角等参数空间对混合的影响规律。而这些将有助于提高人们对混合机理的理解，并有助于推动在上述相关领域的设计创新。因此，本文沿用Avalosse 等[1]和Bertrand 等[2]的模型，将错列角从0 推广到π/2，并建立了相应的物理模型和数学模型。自行开发一种四阶显式Runge-Kutta 积分算法来实现拉格朗日粒子追踪，采用Poincaré 截面、拉格朗日相干结构以及粒子群统计分析等方法来研究分布混合。此外，还比较了不同错列角的功耗。

1 问题描述

如图1所示，存在一个8字形的空腔，放置一对相同的三角形转子，两个转子分别以O1和O2为圆心，逆时针等速旋转。设转子的旋转周期为T，则旋转频率N=1/T。两转子的相对位置如图1所示，空腔壁的外半径为R，转子中心距为C，当O1A1与O2A2平行时，错列角为0。图1对应的是O1A1处于铅锤位置的情况，定义此时刻初始时间。因此，初始时刻的错列角也就等于O2A2与过O2的铅垂线之间的夹角。由于两个转子等速旋转，因此，错列角时刻保持不变。

转子的转速为0.5 r/min，故转动周期T为120 s，频率为N= 1/120 s-1。几何参数如下：流域外径R =30 mm，流域内径R0= 10 mm，中心距C = 48 mm，转子边长L= 40 mm。使用四个典型的错列角度来识别不同的混合情况：θ分别为0、π/6、π/3和π/2。

1.1 建模

研究二维高黏牛顿流体的瞬态等温层流问题。做如下假设：

（1）壁面无滑移；

（2）等温层流；

（3）可忽略的重力和离心效应；

（4）流域全充满。

控制方程如下：

式中，ω 为左右转子的角速度，s-1;如前，转子频率N=1/120 s-1；8字形空腔壁面保持静止。

1.2 数值模拟方法

利用有限元分析软件ANSYS Polyflow 17.0 CFD对上述方程定义的二维瞬时流动进行了数值模拟。8字形流域如图1所示，一对半径为R的相贯圆弧封闭的区域减去一对半径为R0的区域，其中，R0比三角形转子内切圆半径小0.2 mm。采用网格叠加技术（MST）来逼近移动边界。如图2所示，流域和转子分别用四边形网格进行网格划分，然后按照时间步长进行叠加处理。采用ICEM17.0来划分网格。利用对称性和周期性来降低计算成本及提高网格密度，只模拟T/3时间跨度，因此使用上述瞬态模型进行数值模拟，只需要将速度场求解推进到40 s即可。

流动域的左侧以500×80进行网格划分，右侧与左侧相同，因此，在周向和径向上，每个部分分别有1000和160个网格节点。每个三角转子也用曲边三角形镂空，以保证四边形网格易于铺展，减少网格尺寸及提高网格匹配精度，每个转子在周向和径向的网格尺寸均为360×40，如图2 所示。当转子转过2π/3 时，采用100 个时间步，每个时间步等于0.4 s，那么0～40 s 时间跨度内具有101 个时间节点。使用网格叠加技术和作者开发的混合动力学计算代码来完成整体计算。当错列角为0 时，流域（红色）和转子（左蓝右绿）网格叠加结果如图2（a）所示。图2（b）为啮合区域局部放大图。

图2 错列角为0时计算模型在t=0处的网格划分Fig.2 Meshing for the computational model at t=0 for a zero staggered angle

将叠加后的网格导入Polyflow 软件中，进行瞬态速度场求解。迭代参数设置如下：时间步长初始值为T/400，时间上限为0.5T。时间积分采用Galerkin 法，时间按照时间步长向前移动。收敛标准为10-5，每个时间步最大迭代步长定义为30 步。最后，速度结果以CSV 格式导出，读入自行开发混合动力学计算程序，进一步开展混合动力学研究。

图3（a）～（d）为t=0时不同错列角下转子的速度矢量分布。为了便于观察，只使用部分网格点来显示速度矢量分布。图3（e）、（f）为错列角π/3和π/2全网各节点局部放大视图。速度矢量图表明，对于四种错列角，转子周围都有高速带，三角形顶点区都有局部速度峰值。空腔壁面附近存在一条8字形低速带。在转子顶点与空腔壁面附近出现了局部低速涡旋。在t=0时，错列角为π/3时，单个双曲旋涡出现在转子连心线中点附近；而错列角为π/2时，一对双曲旋涡出现在连心线的中垂线中点上下位置附近。双曲点的存在对于混沌混合具有重要意义。

图4为当两个转子旋转部分靠得最近的时刻的速度矢量分布。当错列角为0 时，选择t = 10 s 时刻，对于其他错列角，则选择左右转子顶点首次相遇的时刻，由于错列角的不同，对应的时间分别是15、10 和5 s。如图4（a）、（c）显示，在双曲线的顶点中间出现了椭圆涡，类似于图3（d）。单个双曲线涡点出现在连心线的中垂线上，如图4（b）、（d），类似于图3（a）。

图3 在t=0时刻不同错列角的速度矢量分布/(m/s)Fig.3 The velocity vector distribution of different error angles at time t=0

2 粒子追踪与验证

示踪粒子的运动可以通过求解式(5)进行追踪:

式中，x =(x,y) 为示踪粒子的位置矢量；v(x,y,t)为速度矢量；与式（1）相同，时间t 的单位为s。式（5）用来描述无扩散效应的标记不同颜色的同种物料间的混合。当时间t小于40 s时，通过上述数值模拟方法求解速度矢量。当时间t大于40 s时，利用周期性通过式(b)求得任意时刻的速度场：

式中，MOD（,）为求模运算。

图4 不同错列角速度矢量分布/(m/s)Fig.4 Different staggered angular velocity vector profiles

自编程序，采用标准的四阶Runge-Kutta 格式实现了前锋跟踪计算。对于任意不在网格节点及时间节点上一般情况，需要对时间和空间位置进行插值：一方面，对于任意时刻t，采用二次拉格朗日插值法计算速度值时，采用覆盖时刻t 的k-1、k 和k+1这三个相邻的时间节点。这里的时间步长对应于前面的时间步长，在这个时间步长上，采用网格叠加技术（MST）求解欧拉速度场；另一方面，对于空间9 节点二次插值，将笛卡尔坐标系统（x，y）转化为极坐标（r，θ）来加速计算效率[20-24]。为了实现准确的预测，采用相对较小的时间间隔dt=1.0×10-2s(约为8.33×10-5T)完成粒子跟踪计算。

为了验证网格叠加和粒子跟踪技术的准确性，在错列角为0的情况下，与前人的实验结果[1]进行比较。使用图2 所示的网格划分和前述时间间隔dt，检查了转子旋转2π/3 后的示踪剂混合图像，发现实验结果与数值模拟结果有较好的一致性[21]。

3 结果与讨论

3.1 Poincaré截面

将64个不同颜色的粒子均匀地分布在一个1 mm×1 mm的正方形中，随机滴入流域中，进行前锋追踪计算。每隔一个周期记录所有点的位置，计算100个周期来获得Poincaré截面，如图5所示。典型的“KAM岛”出现在转子每边的中间——左红右蓝。“KAM岛”的面积变化并不明显。例如，对于错列角为0的情况，包含p1=(0.00075,-0.0161)的“KAM岛”，经过一个周期时间T 之后，返回到它们原来的位置，从p2= (0.048, -0.0161)来的点也是如此。一方面,这意味着“KAM岛”区域的流体与转子旋转几乎同步。图5（b）一部分流体只局限于转子的本地搅拌区，局限在左腔或右腔，而彼此不互相混合。另一方面，“KAM岛”的流体从未与外面的流体混合。

3.2 有限时间Lyapunov指数和拉格朗日相干结构

图5 不同错列角的Poincaré截面Fig.5 Poincaré cross sections of different staggered angles

为了确定合适的参数，包括确定初始半径h 和初始半圆周上点的个数N。以错列角π/2 组合为例，取一个特定点p3=(0，-0.024)来试算。采用不同的时间间隔值dt 进行前锋追踪,计算得到在t=T/3(40 s)是FTLE 值,如图6 所示。结果发现，在不同的时间间隔dt和不同的半径h下，得到的FTLE 值几乎相同，当点数N 等于或大于16 时，FTLE 值趋近相同。因此，在实际计算中，N = 16，对应初始圆周围分布的32个点，h=0.01 mm,dt=1.0×10-2s。

图6 p3点的FTLE对不同的时间间隔、初始圆半径和分布点数的敏感性Fig.6 Sensitivity of FTLE of p3 using different values of time interval,initial circle radius and point number

对流动域进行0.5 mm 的网格划分，选择网格线的交点计算FTLE 值。图7 为四种不同错列角在有限时间t= T/3 时的FTLE 场分布。FTLE 分布相似，无明显差异。FTLE 值最大的带从左转子的右下峰到右转子的最高点穿过啮合带。次级大FTLE 值带接转子顶点，转子顶点周围出现局部较大的伴随区域。这些条带代表了FTLE 的脊，揭示了拉格朗日相干结构，表明了不同的排斥和吸引作用的存在。另外，与图5相似，在转子边缘底部发现了局部的低值区，这与KAM 群岛一致。此外,一对FTLE 值较低的区域伴随FTLE 峰脊的上下两侧，分别位于左右两个子流域之内，如图7 所示。将拉格朗日相干结构与Poincaré 截面相结合，揭示了流场中混合强度的不同，高效混合与嵌入式常规规则层流混合的KAM岛并存，三角形转子混合仍有提升的空间。

图7 t=T/3时不同错啮角的FTLE值分布Fig.7 The distribution of FTLE value of different staggered angles at t=T/3

进一步进行粒子群混合示踪分析，在以p4=(0,-0.0185)为中心、边长为3 mm的正方形区域内均布10000个流体粒子，进行前锋追踪运算，将时间推进到2T(240 s)，揭示图7中FTLE分布背后的机理。四个不同错列角的情况下，发现了相似的示踪剂混合图像。示踪剂的拉伸发生在跨越相互啮合区的FTLE脊上，这表明了其排斥性能。也发现FTLE脊带提供了一种方式触发混沌混合，经过此处示踪剂反复拉伸分离为两部分，如图8所示。另一方面，示踪剂也发现局部折叠，这是混沌混合发生触发的第二种机理。可见，四种错列角的混沌触发机理相似。

3.3 混合方差指数表征

根据图5 和图7 的结果，确定了两个典型点。例如，p1位于KAM 岛内，而p4=(0,-0.0185)则位于混沌区域。用一个以p1为中心、边长为1 mm的方形示踪料条s1和一个以p4为中心、边长为3 mm 的方形示踪料条s2来展示四个交错角下不同粒子群的混合演化。每个料条均匀分布着400个粒子。进行粒子跟踪计算，追踪时间达9T(1080 s)。

为了定量评价上述颗粒混合，参考作者之前的工作[21-24]。基于分割单元格子的粒子数的方差计算如下:

因此，无量纲方差I的一个方便的指标可以定义为：

图8 起始于p3的示踪剂经过2 T演化图像Fig.8 Two-period evolution picture of tracer released from p3

显然，随着混合程度的提高，方差减小。即使是全局混沌流也不会产生一个均匀的点分布，而是每个粒子都有均等的机会分布在任何一个单元格子中。这种状态被称为完全空间随机性(CSR)[29-30]。对于M ≫1 和N ≫1，方差指数的相应极限可由泊松概率质量函数得到：

将流体区域排出转子的区域，使用边长为0.25 mm的正方形单元格子来划分网格。当转子旋转时，单元格子的分布随时间更新。对于四种错列角0、π/6、 π/3和π/2，对应的单元格子数约为64215、64129、64068和64205个。在给定时刻t，分别更新粒子群位置和单元格子分布，根据式（10）～式（12）进行统计计算。图9（a）显示了从p1开始的粒子群的方差指数随时间的变化。当粒子被困在一个KAM区域时，如图5所示，准周期波动出现，但没有随时间出现明显的渐近衰减。从图9（b）可以看出，即使经过9个完整的转子旋转周期，这些颗粒仍然局限在KAM区域，难以分散。局部放大的视图捕捉到混合模式的细节为四个交错的角度和规则的粒子仍然保持对齐有序。在图8（b）中，当一个粒子最初释放p4，错列角为π/3时，九全期的准周期性的轨迹记录表明一个周期轨道左和右摄像头单独存在。

当粒子群以p4为中心时，图10（a）绘制了初始粒子群的方差指数随时间的衰减，表明方差指数呈指数下降。差异下降率存在于早期阶段，此时错列角π/3出现下降速度超过其他错列角情况，当时间移向9T 时，接近相同的渐近值，略大于CSR 极限(0.0025)，这意味着有KAM岛阻止了流体在全流域范围内的无序混合，不同情况下的最终混合差异在统计学上并不明显。此外,在t=9 T 混合结果记录在图10（b）、（c）和（d），错列角分别为0、π/3、π/2。

3.4 能耗对比

为了比较不同错列角度下的能耗值，利用旋转轴上一米的长度来计算能耗。因此，式（13）成立：

图9 p1为中心的粒子群经过t=9T(1080 s)演化图像及混合方差指数分析Fig.9 Particle evolutions and mixing variance indices until t=9T(1080 s)initially released from the clusters centered at p1

式中，P(t)为W中的瞬时功耗，为双收缩。Ω代表流域。错列角对功耗的影响如图11（a）所示。由于时间的对称性和周期性，图中只显示了T/3(40 s)。不同的波型几乎以相同的振幅出现。对式(13)，每40 s即进行积分。如图11（b）所示，不同错列角的能耗几乎相同。图11（c）、（d）对比了t=0和T/12(10 s)错列角度为0时空腔内的功耗分布，分别对应图11（a）中的最小和最大能耗。一对高功耗区出现在转子的顶点附近，这看起来像蝴蝶，而近零能耗区域被附加到每个转子边缘的中间部分,表明进一步降低能耗的相关修改转子的几何形状的山峰。在今后的优化设计中，应着重研究转子的尖端结构，比如采用光滑的圆弧结构取代尖点以降低功耗，缩小转子中心距来打破KAM岛区的存在，扩大混沌混合区间，尽可能地提高混合性能。关于这一点，我们将进一步开展研究。

4 结 论

本文建立了一对不同错列角三角形转子搅拌下的空腔混合物理模型和数学模型。采用有限元法求解了高黏度牛顿流体的瞬态流场，采用固定网格法对流动域进行网格划分，采用网格迭加技术考虑了内部转子的存在。开发了一个标准的四阶Runge-Kutta 格式来执行粒子跟踪计算。前人的实验结果验证了本文有限元方法的准确性，两者之间有很好的一致性。错列角增加时，从0 到π/2、π/6发现类似的混合动力学机理，混合指数统计反映了几乎一致的定量化混合指标，而且功耗仍然几乎不变。相同的机制支配了腔内的混合，其中部分混沌混合与规则层流混合的KAM 岛混合并存。部分示踪剂被限制在左、右转子的子区，存在与转子几乎同步旋转的规则层流混合区域。相比之下，混沌混合可以通过两种方式实现。一方面，从拉格朗日相干结构分布来看，FTLE脊提供了将流体示踪剂反复拉伸成两部分的机会，从而创造了一种实现混沌混合的方式。另一方面，当流体示踪剂以逆时针方向绕两转子运动时，流体示踪剂发生局部折叠，导致流体界面发生一系列的拉伸、折叠和重定向，成为出发混沌混合的另外一种机制。猜想当中心距不断缩小或者转子外径不断扩大时，转子顶端会扫略KAM 岛区域，导致KAM 岛破碎，尺度会减小，甚至消除，作者将进一步开展研究。