安 琪 黃志剛 胡淑珍 米國(guó)強(qiáng) 李 賀

(1. 北京工商大學(xué)人工智能學(xué)院，北京 100048；2. 塑料衛(wèi)生與安全質(zhì)量評(píng)價(jià)技術(shù)北京市重點(diǎn)實(shí)驗(yàn)室，北京 100048；3. 中國(guó)農(nóng)業(yè)機(jī)械化科學(xué)研究院，北京 100048)

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

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

適應(yīng)性網(wǎng)格技術(shù)依托網(wǎng)格重疊技術(shù)，根據(jù)單元內(nèi)inside場(chǎng)變量的變化梯度自動(dòng)劃分網(wǎng)格，并通過(guò)改變Maxdiv值來(lái)調(diào)節(jié)網(wǎng)格的密度，該技術(shù)既能節(jié)省選擇合適網(wǎng)格加密區(qū)域的時(shí)間，又可選擇網(wǎng)格的密度以保證計(jì)算精度。Bertrand等[12]曾運(yùn)用適應(yīng)性網(wǎng)格技術(shù)模擬了雙螺桿擠出機(jī)二維模型的旋轉(zhuǎn)過(guò)程，并驗(yàn)證了該方法能夠精確計(jì)算出細(xì)小間隙的剪切速率值，但是由于驗(yàn)證場(chǎng)景的片面以及驗(yàn)證參數(shù)的單一，該技術(shù)在計(jì)算結(jié)果的驗(yàn)證上缺乏準(zhǔn)確性。

文章擬討論適應(yīng)性網(wǎng)格在異向雙螺桿擠出機(jī)中仿真應(yīng)用狀況，并與邊界層網(wǎng)格加密技術(shù)的仿真結(jié)果作對(duì)比，以期證明該方法的實(shí)用性和準(zhǔn)確性。

1 適應(yīng)性網(wǎng)格技術(shù)

1.1 inside場(chǎng)

在Polyflow軟件中，inside是構(gòu)建流體區(qū)域網(wǎng)格運(yùn)動(dòng)部件輪廓的場(chǎng)邊量，能夠描述流體區(qū)域是否被運(yùn)動(dòng)部件重疊。某一節(jié)點(diǎn)的inside值為0說(shuō)明該節(jié)點(diǎn)屬于流體區(qū)域；某一節(jié)點(diǎn)的inside值為1則說(shuō)明該節(jié)點(diǎn)屬于運(yùn)動(dòng)部件區(qū)域。

適應(yīng)性網(wǎng)格技術(shù)基于單元內(nèi)inside場(chǎng)變量的變化梯度能確定該單元是否需要被細(xì)化。若某一單元的變化梯度大于給定值，這種單元在計(jì)算中需要細(xì)化網(wǎng)格，以獲得更加精確的運(yùn)動(dòng)部件輪廓形狀。若某一單元的變化梯度小于給定值，則該單元不需要細(xì)化。

1.2 相關(guān)參數(shù)

1.2.1 網(wǎng)格細(xì)化頻率參數(shù)(Netep) Nstep參數(shù)表示每隔Nstep個(gè)時(shí)間步長(zhǎng)使用一次適應(yīng)性網(wǎng)格技術(shù)。

1.2.2 網(wǎng)格細(xì)化次數(shù)(Maxdiv) Maxdiv定義為使用適應(yīng)性網(wǎng)格技術(shù)的次數(shù)。

1.3 細(xì)化方案

以文中使用的規(guī)則六面體網(wǎng)格為例，當(dāng)細(xì)化次數(shù)Maxdiv為1時(shí)，一個(gè)網(wǎng)格最多可以被分為8個(gè)子網(wǎng)格；當(dāng)Maxdiv為2時(shí)，一個(gè)網(wǎng)格最多可以被分為64個(gè)子網(wǎng)格。

在仿真過(guò)程中，隨著運(yùn)動(dòng)部件的旋轉(zhuǎn)，流道內(nèi)所有單元都會(huì)與運(yùn)動(dòng)部件接觸，因此它們都會(huì)被細(xì)化。當(dāng)運(yùn)動(dòng)部件離開(kāi)了細(xì)化位置，這些細(xì)化位置的細(xì)化次數(shù)將減少，用以節(jié)約計(jì)算時(shí)間和計(jì)算資源[13]。

2 數(shù)值模擬

2.1 幾何模型的建立

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

2.1.2 建立三維模型 根據(jù)圖1繪制異向螺紋元件的端面型線，利用Solidworks軟件建立如圖2的三維模型，模擬雙螺桿轉(zhuǎn)動(dòng)過(guò)程，驗(yàn)證雙螺桿間是否存在物理干涉，并保證建立的幾何模型不存在構(gòu)造重疊。

2.2 數(shù)學(xué)模型的建立

2.2.1 基本假設(shè)[16-17]

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

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

(3) 不計(jì)重力、慣性力等遠(yuǎn)小于黏性力的體積力。

(4) 機(jī)筒壁面的轉(zhuǎn)子表面無(wú)滑移。

(5) 熔體為雷諾系數(shù)較小的層流流動(dòng)。

2.2.2 數(shù)學(xué)控制方程 根據(jù)上述基本假設(shè)條件，基于等溫假設(shè)，不考慮能量方程，給出相應(yīng)流體控制方程[18]：

連續(xù)性方程：

?·ν=0，

(1)

表1 螺紋元件的幾何參數(shù)

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

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

動(dòng)量方程：

?·p+?·τ=0，

(2)

本構(gòu)方程：

(3)

式中：

?——哈密爾頓算子；

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

p——壓力，Pa；

η——表現(xiàn)黏度，Pa·s；

τ——應(yīng)力張量，Pa；

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

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

(4)

式中：

η∞——無(wú)限剪切黏度，Pa·s；

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

λ——carreau模型時(shí)間參數(shù)，s；

n——熔體非牛頓指數(shù)。

2.3 網(wǎng)格劃分

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

2.4 邊界條件和轉(zhuǎn)速的設(shè)定

流道區(qū)域的出入口邊界設(shè)置法向力和切向力為0，說(shuō)明出入口流體可以自由流動(dòng)；左右孔邊界法向速度和切向力為0，說(shuō)明流體不可貫穿邊界，屬于滑移邊界；外壁面的邊界設(shè)置為法向速度和切向速度為0，表示邊界無(wú)滑移[21]。轉(zhuǎn)速設(shè)置為40 r/min。

表2 聚乳酸的物性參數(shù)

表3 網(wǎng)格單元數(shù)的數(shù)據(jù)

3 結(jié)果分析

3.1 流道區(qū)域的邊界網(wǎng)格

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

圖3 網(wǎng)格細(xì)化后的端面網(wǎng)格

3.2 inside場(chǎng)分析

圖4為螺紋元件流道入口邊界inside場(chǎng)的變化云圖。由圖4可知，Maxdiv值為0時(shí)，運(yùn)動(dòng)部件輪廓模糊，嚙合區(qū)內(nèi)螺紋元件inside的范圍兩兩交叉重疊，對(duì)描述嚙合區(qū)內(nèi)流體流動(dòng)情況存在一定干擾；Maxdiv值為1時(shí)，螺紋元件輪廓基本清晰，但嚙合區(qū)內(nèi)inside的范圍依然存在交叉的現(xiàn)象；Maxdiv值為2時(shí)，螺紋元件輪廓更加清晰，嚙合區(qū)內(nèi)兩螺桿inside場(chǎng)范圍獨(dú)立分明。由此可見(jiàn)，運(yùn)用適應(yīng)性網(wǎng)絡(luò)技術(shù)可以提高螺紋元件輪廓的光滑度和清晰度。

3.3 剪切速率場(chǎng)分析

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

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

圖4 網(wǎng)格細(xì)化后的inside云圖

如圖6所示，高剪切速率主要存在于雙螺桿嚙合區(qū)間隙處和螺棱與機(jī)筒間隙處這兩個(gè)區(qū)域內(nèi)，間隙處流體的流動(dòng)速度快，間隙尺寸小，產(chǎn)生了較高的速度梯度，因此剪切速率高[22]。隨著Maxdiv值的增大，剪切速率分布范圍逐漸變小，分布情況逐漸變好，說(shuō)明隨著網(wǎng)格細(xì)化次數(shù)的增加，計(jì)算的結(jié)果越來(lái)越精確。

3.4 黏性場(chǎng)分析

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

圖5 軸向平均剪切速率圖

由圖8可知，螺棱與機(jī)筒間隙處以及靠近螺棱處的流場(chǎng)黏度較小，螺槽中的黏度值較大。邊界層網(wǎng)格的黏度云圖與適應(yīng)性網(wǎng)格(Maxdiv值為0)的差別不大。對(duì)比3種適應(yīng)性網(wǎng)格的云圖，得出：隨著Maxdiv值的增大，螺槽內(nèi)黏度的區(qū)域隨之?dāng)U大，說(shuō)明網(wǎng)格細(xì)化次數(shù)的增多，能夠增加黏度計(jì)算的區(qū)域，顯示黏度分布的更多細(xì)節(jié)。

3.5 壓力場(chǎng)分析

圖9為4種仿真結(jié)果的流場(chǎng)壓力分布云圖，4種結(jié)果的壓力分布普遍不均勻，且沿?cái)D出方向呈增長(zhǎng)的趨勢(shì)。流道下部的壓力高于上部，螺棱與機(jī)筒的間隙處出現(xiàn)了高壓帶，且流道下部的嚙合區(qū)出現(xiàn)局部高壓。

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

圖7 軸向平均黏度折線圖

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

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

3.6 計(jì)算資源消耗的分析

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

(5)

式中：

Q——擠出機(jī)的體積流率，m3/s；

α、β、γ——與螺桿幾何參數(shù)有關(guān)的常數(shù)，m；

KQ——與溫度相關(guān)的常數(shù)，m3；

n——螺桿轉(zhuǎn)速，r/s。

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

圖9 流道壓力云圖

表4 不同的加密網(wǎng)格對(duì)于計(jì)算過(guò)程的影響

4 結(jié)論

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

(2) 在流道網(wǎng)格單元數(shù)相近的情況下，適應(yīng)性網(wǎng)格技術(shù)計(jì)算到的高剪切速率值少于邊界層網(wǎng)格；在黏度場(chǎng)上，適應(yīng)性網(wǎng)格的細(xì)節(jié)更多；在壓力場(chǎng)云圖上，兩者的表現(xiàn)相近。相對(duì)于邊界層網(wǎng)格技術(shù)，適應(yīng)性網(wǎng)格技術(shù)提升了0.40倍的計(jì)算精度，說(shuō)明適應(yīng)性網(wǎng)格技術(shù)能夠獲得更準(zhǔn)確的計(jì)算結(jié)果?？偟膩?lái)看，適應(yīng)性網(wǎng)格技術(shù)在計(jì)算精度上具有優(yōu)勢(shì)。

(3) 研究?jī)H在常規(guī)雙頭螺紋元件上運(yùn)用了適應(yīng)性網(wǎng)格技術(shù)，之后可將該技術(shù)運(yùn)用到具有復(fù)雜構(gòu)造的螺桿元件上，進(jìn)一步體現(xiàn)適應(yīng)性網(wǎng)格自動(dòng)劃分網(wǎng)格的優(yōu)勢(shì)。