適應性網格技術在異向雙螺桿擠出機中的仿真

安 琪 黃志剛 胡淑珍 米國強 李 賀

(1. 北京工商大學人工智能學院，北京 100048；2. 塑料衛生與安全質量評價技術北京市重點實驗室，北京 100048；3. 中國農業機械化科學研究院，北京 100048)

隨著電腦技術和計算方法的飛速發展，電腦和有限元法相結合的模擬仿真方法被廣泛應用在雙螺桿擠出機上[1]。吳小娟等[2]為詳細了解螺旋擠壓工藝對聚合物物料性能的影響，利用網格重疊技術，借助于Polyflow軟件數值模擬單螺桿均化段擠壓過程。王曉瑾等[3]對螺桿頭、過渡段的平直段和錐形段以及定型段采用八點六面體網格，計算時利用網格重疊技術，按實際情況組合后進行仿真。網格重疊技術能夠有效降低流體區域網格劃分的復雜程度，但存在運動部件附近計算結果不準確的缺陷[4]。面對更加復雜的流道時，網格重疊技術更凸顯了計算不準確的缺點，鄧霽蘭[5]在對三螺桿擠出機中流場靜態進行網格劃分時，先用區塊將螺桿與螺棱間隙處、嚙合區同流場的其余部分分開，再單獨加大間隙區、嚙合區的網格密度，最后通過停留時間的實驗驗證，確定了合適的網格劃分原則及網格精度。

目前在使用網格重疊技術的仿真中，多數[6-7]會憑借經驗單獨加大間隙區、嚙合區的網格密度，即邊界層網格加密技術，用以提高計算的精度。但該法不適用于齒形元件[8]、非嚙合波形元件[9]、楔形推力面非常規螺紋元件[10]等非常規、不規律的復雜元件。為了解決螺桿元件復雜表面的邊界網格劃分問題，Hetu等[11]提出了一種浸入邊界方法，但由于該方法對螺桿元件邊界網格的位置非常敏感且準確性取決于單元網格的大小，即該方法需大量的計算資源和計算時間才可保證結果的精確，所以該方法缺乏實用性。

適應性網格技術依托網格重疊技術，根據單元內inside場變量的變化梯度自動劃分網格，并通過改變Maxdiv值來調節網格的密度，該技術既能節省選擇合適網格加密區域的時間，又可選擇網格的密度以保證計算精度。Bertrand等[12]曾運用適應性網格技術模擬了雙螺桿擠出機二維模型的旋轉過程，并驗證了該方法能夠精確計算出細小間隙的剪切速率值，但是由于驗證場景的片面以及驗證參數的單一，該技術在計算結果的驗證上缺乏準確性。

文章擬討論適應性網格在異向雙螺桿擠出機中仿真應用狀況，并與邊界層網格加密技術的仿真結果作對比，以期證明該方法的實用性和準確性。

1 適應性網格技術

1.1 inside場

在Polyflow軟件中，inside是構建流體區域網格運動部件輪廓的場邊量，能夠描述流體區域是否被運動部件重疊。某一節點的inside值為0說明該節點屬于流體區域；某一節點的inside值為1則說明該節點屬于運動部件區域。

適應性網格技術基于單元內inside場變量的變化梯度能確定該單元是否需要被細化。若某一單元的變化梯度大于給定值，這種單元在計算中需要細化網格，以獲得更加精確的運動部件輪廓形狀。若某一單元的變化梯度小于給定值，則該單元不需要細化。

1.2 相關參數

1.2.1 網格細化頻率參數(Netep) Nstep參數表示每隔Nstep個時間步長使用一次適應性網格技術。

1.2.2 網格細化次數(Maxdiv) Maxdiv定義為使用適應性網格技術的次數。

1.3 細化方案

以文中使用的規則六面體網格為例，當細化次數Maxdiv為1時，一個網格最多可以被分為8個子網格；當Maxdiv為2時，一個網格最多可以被分為64個子網格。

在仿真過程中，隨著運動部件的旋轉，流道內所有單元都會與運動部件接觸，因此它們都會被細化。當運動部件離開了細化位置，這些細化位置的細化次數將減少，用以節約計算時間和計算資源[13]。

2 數值模擬

2.1 幾何模型的建立

2.1.1 繪制端面型線 雙螺桿螺紋元件的端面型線本質上決定了擠出機的剪切、塑化、混合、輸送等性能，端面型線的合理設計對提升雙螺桿擠出機的整體性能具有重要意義[14]。在AutoCAD軟件中，根據表1的螺紋元件的幾何參數繪制如圖1的異向雙螺桿螺紋元件的端面型線，并嚴格保證螺旋嚙合時不干涉[15]。

2.1.2 建立三維模型 根據圖1繪制異向螺紋元件的端面型線，利用Solidworks軟件建立如圖2的三維模型，模擬雙螺桿轉動過程，驗證雙螺桿間是否存在物理干涉，并保證建立的幾何模型不存在構造重疊。

2.2 數學模型的建立

2.2.1 基本假設[16-17]

(1) 熔體充滿整個流道且不可壓縮。

(2) 流場中熔體處處有相同溫度。

(3) 不計重力、慣性力等遠小于黏性力的體積力。

(4) 機筒壁面的轉子表面無滑移。

(5) 熔體為雷諾系數較小的層流流動。

2.2.2 數學控制方程 根據上述基本假設條件，基于等溫假設，不考慮能量方程，給出相應流體控制方程[18]：

連續性方程：

?·ν=0，

(1)

表1 螺紋元件的幾何參數

圖1 異向螺紋元件的端面型線

圖2 異向螺紋元件的三維模型

動量方程：

?·p+?·τ=0，

(2)

本構方程：

(3)

式中：

?——哈密爾頓算子；

υ——速度矢量，m/s；

p——壓力，Pa；

η——表現黏度，Pa·s；

τ——應力張量，Pa；

D——形變速率張量，s-1。

研究仿真的是聚乳酸在異向雙螺桿擠出機的等溫流動過程，聚乳酸的物性參數見表2[19]。聚乳酸為非牛頓黏彈流體，式(3)遵循式(4)中的Bird-carreau模型[20]，用于表述聚乳酸的流變特性。

(4)

式中：

η∞——無限剪切黏度，Pa·s；

η0——零剪切黏度，Pa·s；

λ——carreau模型時間參數，s；

n——熔體非牛頓指數。

2.3 網格劃分

將流道和異向螺紋元件的幾何模型導入Gambit分別劃分網格，流道網格為六面體網格，網格數為8 640個。徑向網格分段數為8，軸向網格分段數為15，圓周方向網格數為48，不添加邊界層；螺紋元件網格為四面體網格，尺寸為2 mm。流道網格單元數的數據如表3。

2.4 邊界條件和轉速的設定

流道區域的出入口邊界設置法向力和切向力為0，說明出入口流體可以自由流動；左右孔邊界法向速度和切向力為0，說明流體不可貫穿邊界，屬于滑移邊界；外壁面的邊界設置為法向速度和切向速度為0，表示邊界無滑移[21]。轉速設置為40 r/min。

表2 聚乳酸的物性參數

表3 網格單元數的數據

3 結果分析

3.1 流道區域的邊界網格

圖3為異向雙螺桿螺紋元件流道入口邊界的網格劃分情況。由圖3可知，隨著Maxdiv值的增加，流道入口截面的網格越來越細致；邊界層網格的加密區域集中在嚙合區、螺桿間隙等狹小區域。

圖3 網格細化后的端面網格

3.2 inside場分析

圖4為螺紋元件流道入口邊界inside場的變化云圖。由圖4可知，Maxdiv值為0時，運動部件輪廓模糊，嚙合區內螺紋元件inside的范圍兩兩交叉重疊，對描述嚙合區內流體流動情況存在一定干擾；Maxdiv值為1時，螺紋元件輪廓基本清晰，但嚙合區內inside的范圍依然存在交叉的現象；Maxdiv值為2時，螺紋元件輪廓更加清晰，嚙合區內兩螺桿inside場范圍獨立分明。由此可見，運用適應性網絡技術可以提高螺紋元件輪廓的光滑度和清晰度。

3.3 剪切速率場分析

圖5是流道流場的軸向平均剪切速率圖。分別觀察適應性網格技術仿真結果的變化趨勢可以得到，軸向平均剪切速率曲線隨著z軸距離的變化輕微上揚；由于轉速只有40 r/min，所以曲線的上下波動不明顯。因為網格密度的增加，系統能夠計算到存在于雙螺桿的間隙處和螺桿與機筒的間隙處等小區域的高剪切速率值，所以Maxdiv值為2的曲線處在折線圖最上層， Maxdiv值為0的曲線處在最下層，即Maxidiv值越大，計算到軸向平均剪切速率值越大，曲線就越高。

由于邊界層網格的網格加密區域和網格密度集中在螺桿機筒間隙處，所以該網格能夠計算到更多高剪切速率，如圖6(d)邊界層網格的剪切速率云圖所示，同時比Maxdiv值為1的曲線位置更高，處于Maxdiv值為2的曲線與曲線Maxdiv值為1的曲線之間。

圖4 網格細化后的inside云圖

如圖6所示，高剪切速率主要存在于雙螺桿嚙合區間隙處和螺棱與機筒間隙處這兩個區域內，間隙處流體的流動速度快，間隙尺寸小，產生了較高的速度梯度，因此剪切速率高[22]。隨著Maxdiv值的增大，剪切速率分布范圍逐漸變小，分布情況逐漸變好，說明隨著網格細化次數的增加，計算的結果越來越精確。

3.4 黏性場分析

根據式(4)可以得出，黏度隨著剪切速率的增大而減小[23]。對比圖5和圖7可知，兩者曲線的變化趨勢符合規律。由于圖5和圖7的曲線是選取平均值繪制的，分析的是剪切速率和黏度值兩者軸向的宏觀變化，加之網格局部細化帶來的精度升級，所以Maxdiv值為1、曲線Maxdiv值為2和曲線邊界層網格的黏度數值變化在誤差允許的范圍內[24-25]。

圖5 軸向平均剪切速率圖

由圖8可知，螺棱與機筒間隙處以及靠近螺棱處的流場黏度較小，螺槽中的黏度值較大。邊界層網格的黏度云圖與適應性網格(Maxdiv值為0)的差別不大。對比3種適應性網格的云圖，得出：隨著Maxdiv值的增大，螺槽內黏度的區域隨之擴大，說明網格細化次數的增多，能夠增加黏度計算的區域，顯示黏度分布的更多細節。

3.5 壓力場分析

圖9為4種仿真結果的流場壓力分布云圖，4種結果的壓力分布普遍不均勻，且沿擠出方向呈增長的趨勢。流道下部的壓力高于上部，螺棱與機筒的間隙處出現了高壓帶，且流道下部的嚙合區出現局部高壓。

圖6 流道入口截面的剪切速率云圖

圖7 軸向平均黏度折線圖

圖9(a)顯示的壓力分布呈中心發散式的，壓力隨著距離逐漸減小。圖9(b)顯示的壓力分布呈條狀，由于螺棱與機筒間存在間隙，在瞬態模擬中，壓力集中在螺棱處，隨著兩側尺寸的變化，壓力分布向兩側衰減成條狀，與螺槽內分布的低壓力形成鮮明的對比。圖9(c)顯示的壓力分布：在條狀的基礎上，模糊了高低壓的邊界，壓力進一步衰減，最終呈梯田狀。這是由于網格細化次數的增多，網格的尺寸不斷變小，系統能夠更加精密地計算螺棱螺槽等細小尺寸內的壓力大小，壓力分布的輪廓自然變得更加均勻連續。說明運用適應性網格技術能夠獲得精度更高的計算結果，增加壓力分布的細節，有效提高壓力分布云圖的清晰度。圖9(d)顯示的壓力分布接近于圖9(b)，即邊界層網格和適應性網格的壓力分布相近。

圖8 流道入口截面的黏度云圖

3.6 計算資源消耗的分析

體積流率是研究流場流動的重要參數，數值的大小將直接影響到產品的質量，大的流率可提升生產效率，而過大的流率會嚴重影響產品質量。并根據式(5)可知，體積流率主要受溫度、轉速和螺桿幾何參數的影響[26]。同時，由于設置的溫度、轉速和幾何參數條件都是不變的，出口的體積流率在理論上是定值，故以出口流率量計算結果的精確性。根據表4，以Maxdiv值為0的計算時間、最大內存和出口流率作為參考值，Maxdiv值為1的仿真消耗10.94倍的參考計算時間和6.36倍的參考最大內存，得到1.35倍的計算精度；Maxdiv值為2的仿真結果消耗69.49倍的參考計算時間和39.16倍的參考最大內存，得到1.48倍的參考計算精度。對比Maxdiv值為1和2時的數據，后者比前者1多花費58.55倍的參考計算時間和32.80倍的參考最大內存，提高了0.13倍的參考計算精度，故得出：雖然增大Maxdiv值能加密網格，提升計算精度，但對計算資源的消耗也同樣增加；Maxdiv值為2的仿真雖然消耗大量的計算資源，但是計算精度的提升有限，故Maxdiv值為1的仿真更加適用于常規研究中。在運用適應性網格技術的過程中，應當考慮計算時間和計算機的性能，選擇合適的網格細化次數。

(5)

式中：

Q——擠出機的體積流率，m3/s；

α、β、γ——與螺桿幾何參數有關的常數，m；

KQ——與溫度相關的常數，m3；

n——螺桿轉速，r/s。

邊界層網格的仿真消耗4.79倍的參考計算時間和5.24倍的參考最大內存，得到1.08倍的參考計算精度。與網格單元數相近的適應性網格技術(Maxdiv值為1)對比，適應性網格比邊界層網格多花費6.15倍的參考計算時間和1.12倍的參考最大內存，提高了0.40倍的參考計算精度。雖然適應性網格在自動劃分網格上多花費了計算時間，但是在整個操作上節省了選擇網格加密區域和網格密度的時間，故適應性網格的計算時間在可接受的范圍內。兩者在最大內存上的差距并不大。相對于邊界層網格技術，適應性網格技術提升了0.40倍的計算精度，說明適應性網格技術能夠獲得更準確的計算結果。

圖9 流道壓力云圖

表4 不同的加密網格對于計算過程的影響

4 結論

(1) 隨著Maxdiv值的增大，能夠提高剪切速率場、黏度場和壓力場等流場分布云圖的分辨率，使得圖像更加精細，有助于分析螺桿流道中復雜的流動情況。雖然增大Maxdiv值能加密網格，提升計算精度，但對計算資源的消耗也同樣增加；Maxdiv值為2的仿真雖然消耗大量的計算資源，但是計算精度的提升有限，故Maxdiv值為1的仿真更加適用于常規研究中。在運用適應性網格技術過程中，應當考慮計算時間和計算機的性能，選擇合適的網格細化次數。

(2) 在流道網格單元數相近的情況下，適應性網格技術計算到的高剪切速率值少于邊界層網格；在黏度場上，適應性網格的細節更多；在壓力場云圖上，兩者的表現相近。相對于邊界層網格技術，適應性網格技術提升了0.40倍的計算精度，說明適應性網格技術能夠獲得更準確的計算結果。總的來看，適應性網格技術在計算精度上具有優勢。

(3) 研究僅在常規雙頭螺紋元件上運用了適應性網格技術，之后可將該技術運用到具有復雜構造的螺桿元件上，進一步體現適應性網格自動劃分網格的優勢。