不同錯列角三角形轉(zhuǎn)子作用下的分布混合數(shù)值模擬

徐百平，劉彪，劉堯，譚壽再，杜遙雪，劉春太

（1 五邑大學智能制造學部，廣東江門529020； 2 廣東輕工職業(yè)技術(shù)學院廣東省高分子材料先進加工工程技術(shù)研究中心，廣東廣州510300； 3 鄭州大學橡塑模具國家工程研究中心，河南鄭州450002）

引 言

同向雙螺桿擠出機具有混合性能好、自清潔效果好、脫氣效率高等優(yōu)點，在聚合物加工、反應擠出、食品加工以及制藥等領(lǐng)域獲得了廣泛應用。這類機器內(nèi)，一對形狀相同的內(nèi)部轉(zhuǎn)子，安置在水平腔內(nèi)同向等速旋轉(zhuǎn)，其截面由幾段不同半徑的圓弧組成。為了實現(xiàn)自潔功能，避免運動干涉，兩個轉(zhuǎn)子之間有特定的錯列角。一對三角形轉(zhuǎn)子類似于一對三頭雙螺桿元件和捏合塊，在一定程度上與擠出機的捏合盤具有相同的混合原理。這樣的幾何結(jié)構(gòu)，對于螺紋元件的橫截面內(nèi)混合也有一定的參考意義。然而，當需要更大的工作空間和熱交換面積時，最簡單的方法是縮小轉(zhuǎn)子尺寸以提供更多的空間，但是通常會損失部分自清潔特性。例如，當選擇一對三角形轉(zhuǎn)子時，盡管失去了自清洗功能，但可以提供更多的空間和傳熱表面。而且當轉(zhuǎn)子尺寸縮小后，可以改變錯列角來完成不同的混合過程。但這樣的尺寸放大縮小以及錯列角的選擇對混合室內(nèi)的混合效率影響效果還需要進一步研究。

1997 年，Avalosse 等[1]通過實驗觀察和數(shù)值模擬研究了三角形轉(zhuǎn)子對作用下的混合問題。隨后Bertrand 等[2]提出了另一種虛擬內(nèi)嵌邊界方法，使用拉格朗日乘子法來施加移動邊界條件，進一步進行局部自適應網(wǎng)格細化來精確尖角處邊界條件，從而提高了混合演化細節(jié)預測。當前，計算流體動力學（CFD）取得了巨大的進步，為評估和優(yōu)化反應器設(shè)計提供了一種很有效的方法。有限體積法（FVM）、有限元法（FEM）、光滑質(zhì)點流體動力學法（SPH）等數(shù)值方法現(xiàn)已用于模擬混合裝置中的復雜流動[3-20]。其中，基于有限元的網(wǎng)格疊加技術(shù)（MST）取得了巨大的進步，該方法針對流動域和內(nèi)部運動部件分別進行網(wǎng)格劃分，然后根據(jù)運動部件的運動規(guī)律，在每個時間節(jié)點更新運動部件的運動位置。Kokini團隊[13-16]將有限元數(shù)值模擬方法和疊加網(wǎng)格技術(shù)相結(jié)合，開展了食品加工方面的混合機理研究，取得了豐碩的成果。

針對三角形轉(zhuǎn)子，雖然前人的研究在描述混合細節(jié)和提高數(shù)值精度方面取得了積極進展，但對混合的機理理解仍然不夠清晰。沒有采用動力學研究工具，如統(tǒng)計分析、Lyapunov 指數(shù)及Poincaré 截面等手段，也沒有探索轉(zhuǎn)子大小和幾何形狀、錯列角等參數(shù)空間對混合的影響規(guī)律。而這些將有助于提高人們對混合機理的理解，并有助于推動在上述相關(guān)領(lǐng)域的設(shè)計創(chuàng)新。因此，本文沿用Avalosse 等[1]和Bertrand 等[2]的模型，將錯列角從0 推廣到π/2，并建立了相應的物理模型和數(shù)學模型。自行開發(fā)一種四階顯式Runge-Kutta 積分算法來實現(xiàn)拉格朗日粒子追蹤，采用Poincaré 截面、拉格朗日相干結(jié)構(gòu)以及粒子群統(tǒng)計分析等方法來研究分布混合。此外，還比較了不同錯列角的功耗。

1 問題描述

如圖1所示，存在一個8字形的空腔，放置一對相同的三角形轉(zhuǎn)子，兩個轉(zhuǎn)子分別以O(shè)1和O2為圓心，逆時針等速旋轉(zhuǎn)。設(shè)轉(zhuǎn)子的旋轉(zhuǎn)周期為T，則旋轉(zhuǎn)頻率N=1/T。兩轉(zhuǎn)子的相對位置如圖1所示，空腔壁的外半徑為R，轉(zhuǎn)子中心距為C，當O1A1與O2A2平行時，錯列角為0。圖1對應的是O1A1處于鉛錘位置的情況，定義此時刻初始時間。因此，初始時刻的錯列角也就等于O2A2與過O2的鉛垂線之間的夾角。由于兩個轉(zhuǎn)子等速旋轉(zhuǎn)，因此，錯列角時刻保持不變。

轉(zhuǎn)子的轉(zhuǎn)速為0.5 r/min，故轉(zhuǎn)動周期T為120 s，頻率為N= 1/120 s-1。幾何參數(shù)如下：流域外徑R =30 mm，流域內(nèi)徑R0= 10 mm，中心距C = 48 mm，轉(zhuǎn)子邊長L= 40 mm。使用四個典型的錯列角度來識別不同的混合情況：θ分別為0、π/6、π/3和π/2。

1.1 建模

研究二維高黏牛頓流體的瞬態(tài)等溫層流問題。做如下假設(shè)：

（1）壁面無滑移；

（2）等溫層流；

（3）可忽略的重力和離心效應；

（4）流域全充滿。

控制方程如下：

式中，ω 為左右轉(zhuǎn)子的角速度，s-1;如前，轉(zhuǎn)子頻率N=1/120 s-1；8字形空腔壁面保持靜止。

1.2 數(shù)值模擬方法

利用有限元分析軟件ANSYS Polyflow 17.0 CFD對上述方程定義的二維瞬時流動進行了數(shù)值模擬。8字形流域如圖1所示，一對半徑為R的相貫圓弧封閉的區(qū)域減去一對半徑為R0的區(qū)域，其中，R0比三角形轉(zhuǎn)子內(nèi)切圓半徑小0.2 mm。采用網(wǎng)格疊加技術(shù)（MST）來逼近移動邊界。如圖2所示，流域和轉(zhuǎn)子分別用四邊形網(wǎng)格進行網(wǎng)格劃分，然后按照時間步長進行疊加處理。采用ICEM17.0來劃分網(wǎng)格。利用對稱性和周期性來降低計算成本及提高網(wǎng)格密度，只模擬T/3時間跨度，因此使用上述瞬態(tài)模型進行數(shù)值模擬，只需要將速度場求解推進到40 s即可。

流動域的左側(cè)以500×80進行網(wǎng)格劃分，右側(cè)與左側(cè)相同，因此，在周向和徑向上，每個部分分別有1000和160個網(wǎng)格節(jié)點。每個三角轉(zhuǎn)子也用曲邊三角形鏤空，以保證四邊形網(wǎng)格易于鋪展，減少網(wǎng)格尺寸及提高網(wǎng)格匹配精度，每個轉(zhuǎn)子在周向和徑向的網(wǎng)格尺寸均為360×40，如圖2 所示。當轉(zhuǎn)子轉(zhuǎn)過2π/3 時，采用100 個時間步，每個時間步等于0.4 s，那么0～40 s 時間跨度內(nèi)具有101 個時間節(jié)點。使用網(wǎng)格疊加技術(shù)和作者開發(fā)的混合動力學計算代碼來完成整體計算。當錯列角為0 時，流域（紅色）和轉(zhuǎn)子（左藍右綠）網(wǎng)格疊加結(jié)果如圖2（a）所示。圖2（b）為嚙合區(qū)域局部放大圖。

圖2 錯列角為0時計算模型在t=0處的網(wǎng)格劃分Fig.2 Meshing for the computational model at t=0 for a zero staggered angle

將疊加后的網(wǎng)格導入Polyflow 軟件中，進行瞬態(tài)速度場求解。迭代參數(shù)設(shè)置如下：時間步長初始值為T/400，時間上限為0.5T。時間積分采用Galerkin 法，時間按照時間步長向前移動。收斂標準為10-5，每個時間步最大迭代步長定義為30 步。最后，速度結(jié)果以CSV 格式導出，讀入自行開發(fā)混合動力學計算程序，進一步開展混合動力學研究。

圖3（a）～（d）為t=0時不同錯列角下轉(zhuǎn)子的速度矢量分布。為了便于觀察，只使用部分網(wǎng)格點來顯示速度矢量分布。圖3（e）、（f）為錯列角π/3和π/2全網(wǎng)各節(jié)點局部放大視圖。速度矢量圖表明，對于四種錯列角，轉(zhuǎn)子周圍都有高速帶，三角形頂點區(qū)都有局部速度峰值?？涨槐诿娓浇嬖谝粭l8字形低速帶。在轉(zhuǎn)子頂點與空腔壁面附近出現(xiàn)了局部低速渦旋。在t=0時，錯列角為π/3時，單個雙曲旋渦出現(xiàn)在轉(zhuǎn)子連心線中點附近；而錯列角為π/2時，一對雙曲旋渦出現(xiàn)在連心線的中垂線中點上下位置附近。雙曲點的存在對于混沌混合具有重要意義。

圖4為當兩個轉(zhuǎn)子旋轉(zhuǎn)部分靠得最近的時刻的速度矢量分布。當錯列角為0 時，選擇t = 10 s 時刻，對于其他錯列角，則選擇左右轉(zhuǎn)子頂點首次相遇的時刻，由于錯列角的不同，對應的時間分別是15、10 和5 s。如圖4（a）、（c）顯示，在雙曲線的頂點中間出現(xiàn)了橢圓渦，類似于圖3（d）。單個雙曲線渦點出現(xiàn)在連心線的中垂線上，如圖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時，通過上述數(shù)值模擬方法求解速度矢量。當時間t大于40 s時，利用周期性通過式(b)求得任意時刻的速度場：

式中，MOD（,）為求模運算。

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

自編程序，采用標準的四階Runge-Kutta 格式實現(xiàn)了前鋒跟蹤計算。對于任意不在網(wǎng)格節(jié)點及時間節(jié)點上一般情況，需要對時間和空間位置進行插值：一方面，對于任意時刻t，采用二次拉格朗日插值法計算速度值時，采用覆蓋時刻t 的k-1、k 和k+1這三個相鄰的時間節(jié)點。這里的時間步長對應于前面的時間步長，在這個時間步長上，采用網(wǎng)格疊加技術(shù)（MST）求解歐拉速度場；另一方面，對于空間9 節(jié)點二次插值，將笛卡爾坐標系統(tǒng)（x，y）轉(zhuǎn)化為極坐標（r，θ）來加速計算效率[20-24]。為了實現(xiàn)準確的預測，采用相對較小的時間間隔dt=1.0×10-2s(約為8.33×10-5T)完成粒子跟蹤計算。

為了驗證網(wǎng)格疊加和粒子跟蹤技術(shù)的準確性，在錯列角為0的情況下，與前人的實驗結(jié)果[1]進行比較。使用圖2 所示的網(wǎng)格劃分和前述時間間隔dt，檢查了轉(zhuǎn)子旋轉(zhuǎn)2π/3 后的示蹤劑混合圖像，發(fā)現(xiàn)實驗結(jié)果與數(shù)值模擬結(jié)果有較好的一致性[21]。

3 結(jié)果與討論

3.1 Poincaré截面

將64個不同顏色的粒子均勻地分布在一個1 mm×1 mm的正方形中，隨機滴入流域中，進行前鋒追蹤計算。每隔一個周期記錄所有點的位置，計算100個周期來獲得Poincaré截面，如圖5所示。典型的“KAM島”出現(xiàn)在轉(zhuǎn)子每邊的中間——左紅右藍?！癒AM島”的面積變化并不明顯。例如，對于錯列角為0的情況，包含p1=(0.00075,-0.0161)的“KAM島”，經(jīng)過一個周期時間T 之后，返回到它們原來的位置，從p2= (0.048, -0.0161)來的點也是如此。一方面,這意味著“KAM島”區(qū)域的流體與轉(zhuǎn)子旋轉(zhuǎn)幾乎同步。圖5（b）一部分流體只局限于轉(zhuǎn)子的本地攪拌區(qū)，局限在左腔或右腔，而彼此不互相混合。另一方面，“KAM島”的流體從未與外面的流體混合。

3.2 有限時間Lyapunov指數(shù)和拉格朗日相干結(jié)構(gòu)

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

為了確定合適的參數(shù)，包括確定初始半徑h 和初始半圓周上點的個數(shù)N。以錯列角π/2 組合為例，取一個特定點p3=(0，-0.024)來試算。采用不同的時間間隔值dt 進行前鋒追蹤,計算得到在t=T/3(40 s)是FTLE 值,如圖6 所示。結(jié)果發(fā)現(xiàn)，在不同的時間間隔dt和不同的半徑h下，得到的FTLE 值幾乎相同，當點數(shù)N 等于或大于16 時，F(xiàn)TLE 值趨近相同。因此，在實際計算中，N = 16，對應初始圓周圍分布的32個點，h=0.01 mm,dt=1.0×10-2s。

圖6 p3點的FTLE對不同的時間間隔、初始圓半徑和分布點數(shù)的敏感性Fig.6 Sensitivity of FTLE of p3 using different values of time interval,initial circle radius and point number

對流動域進行0.5 mm 的網(wǎng)格劃分，選擇網(wǎng)格線的交點計算FTLE 值。圖7 為四種不同錯列角在有限時間t= T/3 時的FTLE 場分布。FTLE 分布相似，無明顯差異。FTLE 值最大的帶從左轉(zhuǎn)子的右下峰到右轉(zhuǎn)子的最高點穿過嚙合帶。次級大FTLE 值帶接轉(zhuǎn)子頂點，轉(zhuǎn)子頂點周圍出現(xiàn)局部較大的伴隨區(qū)域。這些條帶代表了FTLE 的脊，揭示了拉格朗日相干結(jié)構(gòu)，表明了不同的排斥和吸引作用的存在。另外，與圖5相似，在轉(zhuǎn)子邊緣底部發(fā)現(xiàn)了局部的低值區(qū)，這與KAM 群島一致。此外,一對FTLE 值較低的區(qū)域伴隨FTLE 峰脊的上下兩側(cè)，分別位于左右兩個子流域之內(nèi)，如圖7 所示。將拉格朗日相干結(jié)構(gòu)與Poincaré 截面相結(jié)合，揭示了流場中混合強度的不同，高效混合與嵌入式常規(guī)規(guī)則層流混合的KAM島并存，三角形轉(zhuǎn)子混合仍有提升的空間。

圖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的正方形區(qū)域內(nèi)均布10000個流體粒子，進行前鋒追蹤運算，將時間推進到2T(240 s)，揭示圖7中FTLE分布背后的機理。四個不同錯列角的情況下，發(fā)現(xiàn)了相似的示蹤劑混合圖像。示蹤劑的拉伸發(fā)生在跨越相互嚙合區(qū)的FTLE脊上，這表明了其排斥性能。也發(fā)現(xiàn)FTLE脊帶提供了一種方式觸發(fā)混沌混合，經(jīng)過此處示蹤劑反復拉伸分離為兩部分，如圖8所示。另一方面，示蹤劑也發(fā)現(xiàn)局部折疊，這是混沌混合發(fā)生觸發(fā)的第二種機理?？梢姡姆N錯列角的混沌觸發(fā)機理相似。

3.3 混合方差指數(shù)表征

根據(jù)圖5 和圖7 的結(jié)果，確定了兩個典型點。例如，p1位于KAM 島內(nèi)，而p4=(0,-0.0185)則位于混沌區(qū)域。用一個以p1為中心、邊長為1 mm的方形示蹤料條s1和一個以p4為中心、邊長為3 mm 的方形示蹤料條s2來展示四個交錯角下不同粒子群的混合演化。每個料條均勻分布著400個粒子。進行粒子跟蹤計算，追蹤時間達9T(1080 s)。

為了定量評價上述顆?；旌?，參考作者之前的工作[21-24]?；诜指顔卧褡拥牧Ｗ訑?shù)的方差計算如下:

因此，無量綱方差I(lǐng)的一個方便的指標可以定義為：

圖8 起始于p3的示蹤劑經(jīng)過2 T演化圖像Fig.8 Two-period evolution picture of tracer released from p3

顯然，隨著混合程度的提高，方差減小。即使是全局混沌流也不會產(chǎn)生一個均勻的點分布，而是每個粒子都有均等的機會分布在任何一個單元格子中。這種狀態(tài)被稱為完全空間隨機性(CSR)[29-30]。對于M ?1 和N ?1，方差指數(shù)的相應極限可由泊松概率質(zhì)量函數(shù)得到：

將流體區(qū)域排出轉(zhuǎn)子的區(qū)域，使用邊長為0.25 mm的正方形單元格子來劃分網(wǎng)格。當轉(zhuǎn)子旋轉(zhuǎn)時，單元格子的分布隨時間更新。對于四種錯列角0、π/6、 π/3和π/2，對應的單元格子數(shù)約為64215、64129、64068和64205個。在給定時刻t，分別更新粒子群位置和單元格子分布，根據(jù)式（10）～式（12）進行統(tǒng)計計算。圖9（a）顯示了從p1開始的粒子群的方差指數(shù)隨時間的變化。當粒子被困在一個KAM區(qū)域時，如圖5所示，準周期波動出現(xiàn)，但沒有隨時間出現(xiàn)明顯的漸近衰減。從圖9（b）可以看出，即使經(jīng)過9個完整的轉(zhuǎn)子旋轉(zhuǎn)周期，這些顆粒仍然局限在KAM區(qū)域，難以分散。局部放大的視圖捕捉到混合模式的細節(jié)為四個交錯的角度和規(guī)則的粒子仍然保持對齊有序。在圖8（b）中，當一個粒子最初釋放p4，錯列角為π/3時，九全期的準周期性的軌跡記錄表明一個周期軌道左和右攝像頭單獨存在。

當粒子群以p4為中心時，圖10（a）繪制了初始粒子群的方差指數(shù)隨時間的衰減，表明方差指數(shù)呈指數(shù)下降。差異下降率存在于早期階段，此時錯列角π/3出現(xiàn)下降速度超過其他錯列角情況，當時間移向9T 時，接近相同的漸近值，略大于CSR 極限(0.0025)，這意味著有KAM島阻止了流體在全流域范圍內(nèi)的無序混合，不同情況下的最終混合差異在統(tǒng)計學上并不明顯。此外,在t=9 T 混合結(jié)果記錄在圖10（b）、（c）和（d），錯列角分別為0、π/3、π/2。

3.4 能耗對比

為了比較不同錯列角度下的能耗值，利用旋轉(zhuǎn)軸上一米的長度來計算能耗。因此，式（13）成立：

圖9 p1為中心的粒子群經(jīng)過t=9T(1080 s)演化圖像及混合方差指數(shù)分析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)。不同的波型幾乎以相同的振幅出現(xiàn)。對式(13)，每40 s即進行積分。如圖11（b）所示，不同錯列角的能耗幾乎相同。圖11（c）、（d）對比了t=0和T/12(10 s)錯列角度為0時空腔內(nèi)的功耗分布，分別對應圖11（a）中的最小和最大能耗。一對高功耗區(qū)出現(xiàn)在轉(zhuǎn)子的頂點附近，這看起來像蝴蝶，而近零能耗區(qū)域被附加到每個轉(zhuǎn)子邊緣的中間部分,表明進一步降低能耗的相關(guān)修改轉(zhuǎn)子的幾何形狀的山峰。在今后的優(yōu)化設(shè)計中，應著重研究轉(zhuǎn)子的尖端結(jié)構(gòu)，比如采用光滑的圓弧結(jié)構(gòu)取代尖點以降低功耗，縮小轉(zhuǎn)子中心距來打破KAM島區(qū)的存在，擴大混沌混合區(qū)間，盡可能地提高混合性能。關(guān)于這一點，我們將進一步開展研究。

4 結(jié) 論

本文建立了一對不同錯列角三角形轉(zhuǎn)子攪拌下的空腔混合物理模型和數(shù)學模型。采用有限元法求解了高黏度牛頓流體的瞬態(tài)流場，采用固定網(wǎng)格法對流動域進行網(wǎng)格劃分，采用網(wǎng)格迭加技術(shù)考慮了內(nèi)部轉(zhuǎn)子的存在。開發(fā)了一個標準的四階Runge-Kutta 格式來執(zhí)行粒子跟蹤計算。前人的實驗結(jié)果驗證了本文有限元方法的準確性，兩者之間有很好的一致性。錯列角增加時，從0 到π/2、π/6發(fā)現(xiàn)類似的混合動力學機理，混合指數(shù)統(tǒng)計反映了幾乎一致的定量化混合指標，而且功耗仍然幾乎不變。相同的機制支配了腔內(nèi)的混合，其中部分混沌混合與規(guī)則層流混合的KAM 島混合并存。部分示蹤劑被限制在左、右轉(zhuǎn)子的子區(qū)，存在與轉(zhuǎn)子幾乎同步旋轉(zhuǎn)的規(guī)則層流混合區(qū)域。相比之下，混沌混合可以通過兩種方式實現(xiàn)。一方面，從拉格朗日相干結(jié)構(gòu)分布來看，F(xiàn)TLE脊提供了將流體示蹤劑反復拉伸成兩部分的機會，從而創(chuàng)造了一種實現(xiàn)混沌混合的方式。另一方面，當流體示蹤劑以逆時針方向繞兩轉(zhuǎn)子運動時，流體示蹤劑發(fā)生局部折疊，導致流體界面發(fā)生一系列的拉伸、折疊和重定向，成為出發(fā)混沌混合的另外一種機制。猜想當中心距不斷縮小或者轉(zhuǎn)子外徑不斷擴大時，轉(zhuǎn)子頂端會掃略KAM 島區(qū)域，導致KAM 島破碎，尺度會減小，甚至消除，作者將進一步開展研究。