復合固體推進劑雙螺桿擠出成型過程DEM-CFD 耦合仿真

王 建，劉 皓，趙亞風，喬曉林，李興剛，趙 慧

（1.北京化工大學機電工程學院，北京 100029；2.北京航天試驗技術(shù)研究所，北京 100074；3.西安航天化學動力有限公司，陜西 西安 710025）

1 引言

雙螺桿擠出成型是火炸藥柔性制造的核心技術(shù)，具有連續(xù)生產(chǎn)進而提高生產(chǎn)效率、減少工藝環(huán)節(jié)進而提高安全程度等優(yōu)勢［1-2］，在固體推進劑生產(chǎn)中的應用將改變傳統(tǒng)生產(chǎn)模式。近年來，數(shù)值模擬技術(shù)在火炸藥螺壓擠出加工生產(chǎn)過程中應用日益突出［3-5］，在安全可靠性方面可發(fā)揮重要作用。但是，由于火炸藥組分的復雜性，相關(guān)模擬技術(shù)［3-5］均把火炸藥物料體系簡化為流體，且只選取了部分較短螺桿結(jié)構(gòu)進行模擬，仿真模擬的精度仍然很差。復合固體推進劑以鋁粉和高氯酸銨氧化劑為主要組分［6-8］，由于固體顆粒含量較多，通過單純的計算流體力學（CFD）的模擬方式難以確保計算精度。復合固體推進劑的過程仿真必須考慮固態(tài)顆粒的運動規(guī)律。離散單元法（DEM）可實現(xiàn)顆粒散狀物料性質(zhì)的有效表征及其加工生產(chǎn)過程的可靠仿真［9-10］。DEM 與CFD耦合可實現(xiàn)復合固體推進劑的過程仿真，并能獲得更高的仿真精度。目前，隨著ANSYS 等計算機輔助分析軟件的發(fā)展，DEM 與CFD 耦合仿真在工業(yè)界得到了廣泛應用［11-14］，但用于火炸藥的雙螺桿加工過程仿真的報道仍然很少。因為復合固體推進劑的雙螺桿加工生產(chǎn)過程仿真不但需要考慮物料顆粒的特性以及流體相的流動特性，還要考慮固液兩相之間的關(guān)系，通過DEM 與CFD 耦合實現(xiàn)復合固體推進劑加工過程的仿真仍然難度較大。

為了推動DEM 與CFD 耦合技術(shù)在固體推進劑雙螺桿加工生產(chǎn)過程仿真中的應用，針對以鋁粉和高氯酸銨氧化劑為主要組分的復合固體推進劑，采用DEM仿真與安息角測試相結(jié)合的標定技術(shù)［15-17］獲得鋁粉、高氯酸銨及其顆?；旌衔锏慕佑|模型參數(shù)，通過流變測試實驗確定液相的黏度模型參數(shù)，采用DEM 與CFD耦合實現(xiàn)復合固體推進劑雙螺桿擠出成型過程的仿真。相關(guān)方法與結(jié)果可為固體推進劑的生產(chǎn)加工過程提供技術(shù)參考，為設(shè)備結(jié)構(gòu)與工藝條件的進一步優(yōu)化和生產(chǎn)工藝的安全可靠性評估提供理論支持。

2 試劑與儀器

復合固體推進劑的固相顆粒物主要由鋁粉和高氯酸銨按質(zhì)量比1∶2 組成；復合固體推進劑液相某比例混合物，西安航天化學動力廠提供。復合固體推進劑由固相顆粒和液相混合物按質(zhì)量比5∶1 混合組成。

黏度計，型號Haake RV1，美國賽默飛世爾科技公司生產(chǎn)。

液體控溫循環(huán)器，型號Haake P1-C35P，美國賽默飛世爾科技公司生產(chǎn)。

3 實驗部分

對固體推進劑中添加的液相混合物進行了黏度測試。設(shè)定溫度窗口為10～70 ℃，測定混合液體的黏度變化。采用Z31 轉(zhuǎn)子，Z43 轉(zhuǎn)筒，裝入液體50.5 mL，為轉(zhuǎn)筒體積2/3。降溫及加熱過程通過液體控溫控制器控制，升溫速率2 ℃·min-1，測溫點保溫30 min，溫度精度為±0.1 ℃，定剪切速率0.1 s-1。圖1 顯示了復合固體推進劑液相混合物的黏度隨溫度的變化曲線。隨著溫度的升高，物料的黏度下降。在仿真計算時，假設(shè)液相為牛頓流體。結(jié)合實際生產(chǎn)工藝，取溫度40 ℃時液相混合物對應的黏度值1.78 Pa·s 用于DEM-CFD耦合仿真。

圖1 混合液體的黏度變化Fig.1 Mixed liquid viscosity changes with temperature

4 模擬仿真

4.1 DEM 仿真

（1）幾何模型

仿真計算的幾何模型為雙螺桿裝配模型，包含機筒、兩根螺桿以及進料口。機筒一端為開口狀態(tài)，進料口和出料口（機筒開口）為常壓口。雙螺桿中心距43 mm，螺桿外徑為52 mm，由螺紋元件、嚙合元件組成，依次包括三塊螺距64 mm 的螺紋元件、一塊螺距32 mm 的螺紋元件、三組60°錯列角厚32 mm 的嚙合元件、五塊螺距32 mm 的螺紋元件、一塊螺距32 mm的反向螺紋元件，總長512 mm，如圖2a 所示。螺桿之間的間隙以及外螺紋與機筒之間的間隙均為1.2 mm。補充幾何平面約束來完善DEM 中模型的處理，如圖2b 所示，添加標紅的平面作為密封面，添加標紅無底長方體盒子作為顆粒工廠。根據(jù)工藝參數(shù)定義螺桿幾何體轉(zhuǎn)動軸及轉(zhuǎn)動方向與轉(zhuǎn)速。研究中涉及兩根螺桿的同向轉(zhuǎn)動，過大的網(wǎng)格數(shù)量使計算機負荷增大，最終確定用于計算的螺桿網(wǎng)格尺寸為2 mm，總網(wǎng)格數(shù)721548。圖2c 顯示了螺桿網(wǎng)格劃分情況。

圖2 雙螺桿擠出機構(gòu)模型Fig.2 Model of the twin-screw extruder

（2）物料模型

鋁粉顆粒泊松比為0.3，密度2700 kg·m-3，剪切模量2.6×1010Pa；高氯酸銨顆粒泊松比0.4，密度1500 kg·m-3，剪切模量1×108Pa。物料模型中包含多種顆粒模型，離散元計算的基礎(chǔ)是顆粒的運動。需要對顆粒形狀和大小等進行定義。將顆粒定義為球形顆粒，標準顆粒半徑為0.1 mm，考慮到真實的鋁粉和高氯酸銨顆粒直徑分布不均勻，設(shè)置顆粒工廠中粒徑在一定范圍內(nèi)隨機產(chǎn)生，最大粒徑為0.3 mm，最小粒徑為0.02 mm，在顆粒工廠中的隨機位置產(chǎn)生。物料最適宜的表面能為5 J·m-2，相關(guān)用于DEM 分析的接觸模型參數(shù)如表1 所示［18］。

表1 接觸模型參數(shù)Table 1 Parameters of the contact model

（3）DEM 仿真控制方程

考慮到研究對象為高固體含量的推進劑，認為材料混合后固體顆粒被包裹在膠體之中，顆粒與顆粒之間的作用力不僅存在于法向作用力，也有較為明顯的切向作用力。Hertz-Mindlin with JKR 模型［20］優(yōu)勢在于存在法向力以及切向粘結(jié)力，無論對于前期的粉體料混合過程還是后期加入凝膠后的漿體推進劑，Hertz-Mindlin with JKR 模型均可較好的解釋實際固體推進劑的物理行為。因此，選擇DEM 中Hertz-Mindlin with JKR 模型作為顆粒接觸力模型。模型的法向作用力取決于顆粒之間的重疊量和表面能，且存在式（1）關(guān)系：

式中，F(xiàn)JKR為法向作用力，N；E為楊氏模量，Pa；R為當量半徑，mm；a為接觸半徑，mm，此處與物理半徑保持一致；σ為顆粒重疊量，mm；η為表面能，J·m-2。

Hertz-Mindlin with JKR 模型提供了吸引凝聚力，即使顆粒并不是直接接觸，顆粒間仍然存在非零凝聚力的最大間隙為σ0，可以通過式（3）計算：

式中，σ0為存在凝聚力的最大間隙，mm；當σ0＜σ時，模型計算返回值為0，此時物料顆粒之間距離過大，不存在凝聚力。

（4）工藝參數(shù)設(shè)置

根據(jù)實際生產(chǎn)情況，螺桿轉(zhuǎn)速統(tǒng)一設(shè)置為70 r·min-1，兩根螺桿同向轉(zhuǎn)動。在DEM 固體顆粒輸送仿真過程中，鋁粉與高氯酸銨質(zhì)量比為1∶2 進料，顆粒按質(zhì)量流率生成。設(shè)置了5 組不同喂料量的實驗，如表2所示。

表2 顆粒喂料量設(shè)定Table 2 Feeding rate of the powders kg·h-1

（5）求解器設(shè)置

考慮后續(xù)與CFD 軟件耦合時時間步長之間的聯(lián)系，將設(shè)定時間步長為2.5e-07 s，計算得到瑞利時間步長的22.9%。經(jīng)過多次試算，在螺桿轉(zhuǎn)速70 r·min-1時，經(jīng)過15 s，顆?？梢詮挠覀?cè)出料口擠出，再經(jīng)過5 s之后發(fā)現(xiàn)顆粒的運動分布情況以及基本趨于穩(wěn)定，最終設(shè)定總時長20 s 進行仿真計算。

4.2 CFD-DEM 仿真

（1）模型與網(wǎng)格劃分

DEM-CFD 耦合仿真中對筒體內(nèi)壁面、進料口、出料口、螺桿區(qū)域進行面劃分，對流道區(qū)域采用填充的體網(wǎng)格劃分，總網(wǎng)格數(shù)483002，網(wǎng)格模型如圖3 所示。因為計算工作量問題，CFD 模擬采用的螺桿模型只選取了螺紋輸送段和嚙合段，省去了DEM 模擬中嚙合段之后的輸送段。采用動網(wǎng)格Porfile 文件指定網(wǎng)格運動。由工藝條件設(shè)定螺桿轉(zhuǎn)動速度與DEM 設(shè)置相同，為70 r·min-1。利用Profile 文件中的表示角速度的語句定義，采用瞬態(tài)方式，取時間0、0.1、1000 s 三個數(shù)據(jù)點，三個數(shù)據(jù)點沿x軸轉(zhuǎn)速均為7.33 rad·s-1。代表其勻速運動，無轉(zhuǎn)動加速度。

圖3 DEM-CFD 仿真螺桿網(wǎng)格劃分Fig.3 Screw mesh generation in DEM-CFD

（2）DEM-CFD 耦合仿真控制方程

顆粒流與流體進行耦合計算的主要模型有兩種：拉格朗日（Lagranginan）模型和歐拉（Eulerian）模型。兩種模型的選擇主要取決于固體含量的多少。拉格朗日模型中只包括固液兩相之間的動量交換，通常在固含量低于10%時選用該模型。對于本研究中復合固體推進劑，固含量為85%，故采用歐拉模型。歐拉模型考慮了固液間的動量交換以及固體顆粒對于流體流動狀態(tài)的影響。歐拉方程是拉格朗日方程的修正方程，在描述耦合模型時在控制方程中添加流體體積分數(shù)φ，當φ=1 時，模型為拉格朗日模型，物料中無固相存在，當φ＜1 時，模型為歐拉模型。耦合后，流體相滿足連續(xù)性方程及動量守恒方程，顆粒系統(tǒng)則遵循牛頓第二定律。

流體連續(xù)性方程為：

動量守恒方程為：

顆粒運動方程滿足：

式中，ρl為液體密度，kg·m-3；φ為液相體積分數(shù)；vl為液相速度，m·s-1；g 為重力加速，9.81 m·s-2；p為壓力，Pa；μ為液體黏度，Pa·s；S為兩相動量交換源項，kg·m·s-1；V為CFD 網(wǎng)格單元體積，m3；Fd，i表示作用在網(wǎng)格單元的流體阻力，N；∑F為顆粒所受的合力，N。

（3）耦合計算條件

首先設(shè)置DEM 模型以及物料參數(shù)，通過讀取UDF 耦合文件進行耦合。DEM 與CFD 軟件進行耦合后，CFD 軟件自動讀取DEM 中的各項參數(shù)作為計算參數(shù)。計算方式為瞬態(tài)計算，重力設(shè)定為z軸9.81 m·s-2，與DEM 中設(shè)定保持一致。螺桿轉(zhuǎn)速為70 r·min-1，流體進入流量設(shè)定為15 kg·h-1，固體顆粒進料質(zhì)量流量設(shè)置為75 kg·h-1。出口設(shè)置為常壓出口。溫度設(shè)定為40 ℃，流體粘度為常數(shù)黏度1.78 Pa·s，密度為920 kg·m-3。

求解器設(shè)置主要注意時間步長。耦合計算需要根據(jù)DEM 中的求解器關(guān)聯(lián)設(shè)置，CFD 耦合計算時的時間步長為DEM 計算時間步長的100 倍，最大迭代次數(shù)為40 次，總時間步10000 次。

5 結(jié)果與討論

5.1 DEM 模擬結(jié)果分析

（1）總進料量選擇

在DEM 模擬中，設(shè)置五組混合進料來判斷顆粒工廠產(chǎn)生顆粒的適宜的喂料速率，分別為30、45、60、75、90 kg·h-1，螺桿轉(zhuǎn)速均為70 r·min-1，其它顆粒屬性設(shè)置與計算邊界條件設(shè)定都相同，物料填充結(jié)果如圖4 所示。圖4a 中，物料主要分布在螺桿的機筒壁面以及螺桿嚙合處，綜合圖4a 中的粒子分布，不難發(fā)現(xiàn)此時進料量過低。分析圖4 中物料分布，低喂料量進料情況下，固體顆粒在進料過程較為分散，表現(xiàn)為小顆粒的局部團聚，物料整體不連續(xù)。這是由于物料進料量過低，顆粒較少，多數(shù)顆粒與顆粒之間距離較大，遠遠大于JKR 計算模型中存在粒子間粘結(jié)力的最大距離［20］。另一現(xiàn)象為物料容易沉降在機筒底部，運動緩慢。這是由于物料受到重力作用的影響，與周圍顆粒粘結(jié)后沉降在螺桿底部，而顆粒的堆積高度小于螺桿與機筒之間的間隙，螺桿的旋轉(zhuǎn)產(chǎn)生的對物料的推動力、摩擦力、剪切力等無法作用到團聚的顆粒上。只有在堆積高度大于螺桿與機筒間隙時才會有一小部分向前運動。這也是靠近進料口的物料明顯多于靠近出料口的螺桿段的物料的原因。圖4b 與圖4c 中，物料隨著進料的增加在螺桿上分布有了較大改善。圖4d 與圖4e中進料口都存在物料堆積，滿足了顆粒之間相互接觸、連續(xù)進料的要求。鑒于過大的進料會給設(shè)備造成較大的負荷，因此最終選用的適宜固體顆粒進料量為75 kg·h-1。

圖4 不同顆粒喂料量條件下在4 s 時刻的物料顆粒分布Fig.4 Material powder distribution at 4 s for different feeding rate of powders

（2）顆粒分布

為確保固體推進劑雙螺桿擠出成型過程仿真在DEM 中的正確性，將仿真結(jié)果進行處理。分析了達到穩(wěn)定狀態(tài)時候，物料在雙螺桿中的分布狀況如圖5a 所示。耿孝正［19］在1993 年研制成功國內(nèi)第一臺可視化雙螺桿擠出實驗裝置，通過此裝置做了固體粒料的雙螺桿輸送實驗，圖5b 為螺桿粒料輸送的實際分布情況。觀察仿真結(jié)果圖5a，紅框區(qū)域分別與圖5b 中視角對應，可以發(fā)現(xiàn)DEM 仿真計算的結(jié)果與耿孝正聚丙烯顆粒料輸送結(jié)果描述的分布情況相符。在上輸送區(qū)內(nèi)物料聚集，在左輸送區(qū)內(nèi)物料堆積偏左下方堆積。由于螺桿的上旋，旋轉(zhuǎn)的螺桿對物料的摩擦拖拽，使得物料有上爬行為。當重力大于摩擦力時，物料下滑，最終在斜面位置保持平衡。螺桿的正右側(cè)方向基本上看不到物料，右輸送區(qū)的物料運動在下嚙合區(qū)處，在機筒隆起的部受阻。DEM 仿真結(jié)果與可視化雙螺桿固體混合試驗結(jié)果吻合，證明了DEM 仿真結(jié)果的正確。因為傳統(tǒng)計算流體力學軟件等對雙螺桿的模擬無法得出固體分布情況。對比這一情況，反映出利用DEM 方法分析雙螺桿固體粒料的輸送是可行的，也證明了DEM中內(nèi)嵌的物理模型適用性較強，在參數(shù)設(shè)置得當?shù)那闆r下，可以計算出接近雙螺桿試驗真實情況的結(jié)果。

圖5 固體顆粒物料在雙螺桿中的分布Fig.5 Solid powder material distribution in the twin-screw extruder

（3）固體粒料在雙螺桿中不同位置的填充率

通過總時長20 s 的計算，在螺桿的徑向方向劃分150 個統(tǒng)計網(wǎng)格單元組，統(tǒng)計每個網(wǎng)格中的粒子數(shù)目。結(jié)合螺桿的位置關(guān)系，統(tǒng)計了填充率隨時間的變化關(guān)系，如圖6a 所示。特別分析了16，17，18，19 s 的雙螺桿的固體顆粒的填充情況，如圖6b 所示。通過分析填充率分布圖，可以發(fā)現(xiàn)，隨著時間推移，嚙合區(qū)填充率逐漸升高。在19 s 時，嚙合區(qū)和螺桿末端反向螺紋元件位置處的物料填充率接近100%；在輸送端位置處的物料填充率較低，大概為20%左右；而在嚙合區(qū)后的螺紋輸送區(qū)域的物料填充率較大，約為30%。隨著時間的推移，螺桿前半部分螺紋輸送端的物料填充率變化不大，而經(jīng)過嚙合段后的螺紋輸送區(qū)域的物料填充率逐漸上升。這是由于在到達16 s 之前，前部分螺桿螺紋輸送段的物料分布已經(jīng)處于穩(wěn)定狀態(tài)；在時間達到19 s 后，螺桿嚙合去后的螺紋輸送段的物料填充率也逐漸達到了穩(wěn)定狀態(tài)。機筒各處填充率不同是由于不同的螺桿構(gòu)型對顆粒的作用力不同。不同于螺桿螺紋輸送段螺紋元件的摩擦拖拽，固體顆粒在嚙合區(qū)的輸送動力主要來源于嚙合塊的擠壓。在嚙合段的物料只有一個沿螺桿軸向的初速度，不存在直接向前運輸?shù)膭恿?，受到嚙合塊的擠壓作用，容易在嚙合區(qū)處形成堆積。而反向螺紋元件位置處的物料顆粒速度受到螺桿作用改變運動方向，也會產(chǎn)生大量物料堆積。

圖6 沿雙螺桿的物料填充率隨時間變化情況Fig.6 Filling rate along the twin-screw at different time

5.2 DEM-CFD 耦合模擬結(jié)果分析

（1）填充率分布

經(jīng)瞬態(tài)分析發(fā)現(xiàn)19 s 物料在螺桿中的輸送已達到穩(wěn)定狀態(tài)，DEM-CFD 耦合模擬均采用19 s 時刻DEM 穩(wěn)定狀態(tài)的結(jié)果進行。圖7 為DEM 法和DEM-CFD 耦合法模擬得到的固體顆粒在螺桿中的填充率分布曲線對比。在相同的工藝條件下，與DEM 得到的結(jié)果相比較，發(fā)現(xiàn)加入質(zhì)量分數(shù)為15%的液相成分后，固體顆粒的分布填充情況大為改善，填充速率更快。前段螺紋元件輸送段的填充率由20% 提升到40%，提升了將近一倍。結(jié)果反映液相的加入極大地改善了固體推進劑的流動性。嚙合段的填充率接近為100%，加入液相后嚙合段向下一段的輸送能力明顯得到了提高。

圖7 沿螺桿方向的填充率分布圖對比Fig.7 Comparison of filling rate distribution diagram in screw direction

（2）顆粒平均速度

圖8a 為顆粒平均速度沿螺桿方向的分布曲線。加入液相的固體推進劑顆粒的平均速度遠大于純固體顆粒的平均速度。整體顆粒的平均速度由DEM 模擬得到的0.1 m·s-1提升到DEM-CFD 耦合得到的0.25 m·s-1。未加液相的固體推進劑顆粒與螺桿的摩擦主要為滑動摩擦，螺桿較難帶動物料運輸；顆粒的輸送主要依靠螺桿的推力。而由液相包裹的顆粒形成漿體，有利于物料對螺桿的黏附，受到的螺桿剪切和摩擦力相對變大，利于物料的輸送，運動速度也相對變大。而且，物料在嚙合處受到嚙合塊的擠壓作用，速度進一步上升。在進料口位置處，耦合計算時的物料平均速度突然上升，可能是由于進料量過大，在進料口大小不變的情況下，液體要達到設(shè)定進料速度必須有足夠大入口速度。圖8b 顯示了顆粒最大速度沿螺桿方向的分布曲線。顆粒最大速度分布規(guī)律并不明顯，耦合前后速度峰值變化不大，但是DEM-CFD 耦合計算得到的最大速度的平均值高于耦合前DEM 計算得到的最大速度平均值，說明液相的加入對顆粒最大速度的提高也有促進作用。

圖8 沿螺桿方向的顆粒速度分布圖Fig.8 Particle velocity distribution along the screw direction

（3）螺桿受力情況

圖9 為DEM 和DEM-CFD 耦合計算得到的螺桿受力分布云圖。圖10 為DEM 和DEM-CFD 耦合計算得到的兩根螺桿分別沿螺桿方向的受力曲線。兩螺桿受力情況基本一致，最大受力點都出現(xiàn)于嚙合區(qū)域及嚙合塊與螺紋元件連接位置處。加入液相前后螺桿受力變化不大。雖然液相的加入增加了物料的流動性，理論上能夠降低螺桿的受力；但是較好的流動性也使得物料在螺桿中的填充率變大，物料堆積更多。因此，流動性的增強和物料的增多對螺桿產(chǎn)生的影響幾乎相互抵消。

圖10 沿螺桿方向的螺桿受力曲線Fig.10 Screw force curves along screw direction

6 結(jié)論

為了解決復合固體推進劑雙螺桿擠出仿真中固液相性質(zhì)差距較大的問題，首先通過DEM 方法對固體推進劑的主要固相顆粒物（鋁粉和高氯酸銨）在雙螺桿中的擠出過程進行了仿真分析，然后用DEM-CFD 耦合的方法對固液兩相流的復合固體推進劑在雙螺桿中的擠出過程進行了仿真分析，與DEM 計算得到的純固體顆粒擠出過程進行對比，得到以下結(jié)論：

（1）對復合固體推進劑固相顆粒在雙螺桿中擠出過程的DEM 模擬分析發(fā)現(xiàn)，DEM 法對于固體推進劑粉料的擠出混合過程得到的固體顆粒分布規(guī)律與雙螺桿可視化實驗結(jié)果相符，這是傳統(tǒng)有限元分析方法難以實現(xiàn)的?？紤]液相的加入，采用DEM 與CFD 相結(jié)合的耦合計算方法可實現(xiàn)更準確的模擬。

（2）DEM 模擬分析了不同喂料量條件下固相顆粒在螺桿中的分布情況，確定了螺桿轉(zhuǎn)速70 r·min-1條件下最優(yōu)的固體顆粒喂料量為75 kg·h-1。

（3）采用流變測試得到40 ℃時復合固體推進劑液相的黏度為1.78 Pa·s，該黏度值用于CFD 及DEM-CFD耦合仿真。在固體顆粒喂料量75 kg·h-1和轉(zhuǎn)速70 r·min-1的工藝條件下，采用DEM-CFD 耦合模擬分析計算了固含量質(zhì)量分數(shù)75%的復合固體推進劑在雙螺桿擠出過程中的填充率、速度和螺桿受力情況。

（4）比對DEM 和DEM-CFD 耦合計算結(jié)果發(fā)現(xiàn)：加入液相物料后，顆粒平均速度由0.1 m·s-1提升到0.25 m·s-1；液相的加入增加了物料在輸送段的填充率，由20%提升到40%，物料的輸送效率也得到提升；螺桿的受力變化不大；在實際生產(chǎn)中應當著重注意嚙合段螺桿及物料的受力情況。