基于光滑粒子流體動力學（SPH）方法的環(huán)模制粒成型數(shù)值模擬方法研究

■孫萬峰 王 禹 孫 宇 武 凱

(南京理工大學機械工程學院，江蘇南京 210094)

環(huán)模制粒工藝及裝備以其穩(wěn)定的生產能力被廣 泛應用于生物質能源和飼料工業(yè)等領域，目前存在成型能耗高、設備磨損快、制品品質不穩(wěn)定的問題[1]。為了突破這些技術難題，工程上往往采用試驗探索的方法。此外，一旦物料或設備變化，則需要重新試驗，效率低、成本高、效果差。而通過數(shù)值模擬的方法來代替試驗探索則能有效解決這一問題[2]。

目前，成型過程數(shù)值模擬相關研究主要結合有限元法和計算流體力學法進行，忽略松散物料顆粒間的相互作用[3]。然而，環(huán)模制粒致密成型過程實際上是物料在模輥壓力作用下從松散狀態(tài)變?yōu)橹旅軤顟B(tài)的過程，原料顆粒間的相互作用對成型工藝有至關重要的影響[4-5]。此外，環(huán)模制粒過程具有非連續(xù)循環(huán)壓制特征，物料被少量多次壓入模孔中，并被逐漸擠出。因此，必須考慮顆粒間的相互作用和非連續(xù)多次壓制條件，從而使仿真與工程實際更吻合。

單個?？字形锪系某尚瓦^程是能反映環(huán)模制粒工藝特點的最小核心單元。目前對于環(huán)模制粒成型機理的研究，多采用將制粒過程簡化為單孔擠壓等方式進行[6]，試驗既能反映物料在?？變仁軌翰D出的核心過程[7]，又能降低試驗（仿真）成本和難度。

光滑粒子流體動力學（smoothed particle hydrodynamics，SPH）是近20多年來逐步發(fā)展起來的一種無網(wǎng)格方法，其主要優(yōu)點是可以處理基于網(wǎng)格的更大的局部變形[8]。這種方法于1997年引入[9]，最初用于天體物理學中的問題，之后進一步發(fā)展將SPH擴展到固體力學[10]。

文章結合LS-DYNA軟件利用SPH法對單個?？字猩镔|材料非連續(xù)多次壓制過程進行數(shù)值模擬，研究成型過程中壓制力的變化情況和物料的流動規(guī)律，并將仿真結果與已有研究對比，以期為生物質環(huán)模制粒成型過程的數(shù)值模擬提供新的方法。

1 原理分析與模型構建

1.1 環(huán)模制粒成型原理

環(huán)模在電機驅動下帶動物料，進而帶動壓輥同向旋轉[11]，如圖1所示。一般認為制粒室分為3個部分，分別是供料區(qū)、變形壓緊區(qū)和擠壓成型區(qū)[12]。物料在環(huán)模和壓輥的作用下，從供料區(qū)經(jīng)過變形壓緊區(qū)，最后進入擠壓成型區(qū)，進而從?？妆粩D出，制成致密顆粒。

圖1 環(huán)模制粒成型過程及其放大圖

當環(huán)模制粒致密成型進入穩(wěn)定擠壓狀態(tài)時，?？變葔褐茟δＰ腿绻剑?）[13]。

式中：v——物料被壓縮速度（m/s）；

vr——壓輥轉速（m/s）；

d——模輥中心距（mm）；

r——壓輥半徑（mm）；

R——環(huán)模半徑（mm）；

h——物料距環(huán)模工作面的距離（mm）；

γ——物料與壓輥圓心連線與水平線之間的夾角。

1.2 SPH方法

SPH 的核心是插值，所有宏觀變量（如密度、壓力、速度等）都通過離散點處的一組值方便地表示為積分插值計算。數(shù)值域中點x處的場變量f（x）由公式（3）計算相鄰粒子的影響而獲得[8]。

N——影響粒子的數(shù)量。

上述插值積分的倒數(shù)是通過普通推導獲得的，不使用有限差分，沒有網(wǎng)格，這是SPH 相對于常規(guī)數(shù)值方法（有限差分或有限元）的主要優(yōu)勢。通過在每個粒子處分配其自己的平滑長度并允許其隨時間變化，SPH 模擬可根據(jù)當?shù)貤l件自動進行調整。由于粒子在計算過程中可以相互作用然后分離，因此SPH能夠處理非常大的變形[10]。

單孔擠壓過程中，物料從橫截面大的擠壓區(qū)進入橫截面小的?？祝l(fā)生較大的變形，有限元方法會出現(xiàn)網(wǎng)格畸變、計算不收斂等問題，為避免這些問題，采用無網(wǎng)格的SPH方法。

1.3 單孔壓制模型的構建

為了研究物料的受力情況，且減少計算量與難度，對環(huán)模進行簡化，提取出單孔模型[7]。環(huán)模制粒機的關鍵零部件為環(huán)模，其工作面如圖2所示。環(huán)模的工作面是由模孔排列形成的，?？着c?？字g存在空隙。模孔之間空隙的大小由?？椎臋M截面及倒角和環(huán)模的剛度需求所決定的。環(huán)模的開孔率被定義為環(huán)模所有?？椎拿娣e與整個環(huán)模工作面總面積之間的比值。環(huán)模的開孔率通常被設計在30%～40%[7]。

圖2 不同?？着帕蟹绞降沫h(huán)模工作面（部分）

由于?？字g存在空隙，物料在擠壓成型區(qū)由一塊致密的塊體被分成若干個小塊分別被擠入模孔內。不考慮環(huán)模邊緣處的模孔，其他?？兹鐖D2 所示，每個?？字車? 個或6 個對稱分布的模孔，因此，當物料在?？咨戏绞艿綌D壓時，在模孔徑向上，物料的受力是對稱的。因而可以沿著圖2 中的虛線提取出單個?？走M行研究物料的受力情況，同時將虛線內的區(qū)域稱為單孔擠壓區(qū)。

圖2中，X表示?？字行木?，為了計算單孔擠壓區(qū)的大小，借鑒環(huán)模開孔率的定義，將單孔擠壓區(qū)的面積定義公式（5）所示。

式中：S單——單孔擠壓區(qū)的面積（mm2）；

r——?？椎陌霃剑╩m）；

ε——環(huán)模開孔率。

2 數(shù)值模擬

采用動力學數(shù)值模擬軟件LS-DYNA實現(xiàn)的SPH方法，因為它能夠模擬大變形問題，同時采用重啟動或者Dynain文件法可以實現(xiàn)多次擠壓過程，可以更好地模擬真實的單孔擠壓過程。

2.1 實體模型的建立與網(wǎng)格劃分

根據(jù)對單孔擠壓區(qū)的分析，其在徑向上是對稱受力的，因而本研究中將多邊形的單孔擠壓區(qū)簡化為圓形。在數(shù)值模擬中，將模型幾何體設置為軸對稱的3D模型，將?？壮叽缭O置成反芻料常用尺寸，具體尺寸如圖3所示。

圖3 數(shù)值模擬的模型及尺寸

在單孔擠壓過程中，料筒和?？變H有其內表面參與物料的壓制，因此，將料筒和?？椎葞缀误w設置成殼單元，從而在不影響計算精度的情況下，減少計算量和計算成本。

物料使用SPH粒子進行計算，物料在擠壓過程中會發(fā)生大變形，使用SPH法可以避免常規(guī)的有限元法出現(xiàn)網(wǎng)格畸變等問題。SPH粒子生成的方式很多，文章使用LS-DYNA的前后處理軟件LS-PrePost 將已有的Solid單元轉換為SPH粒子，具體如圖4所示。

圖4 網(wǎng)格劃分及生成SPH粒子

粒子與?？字g的接觸設置為自動點-面接觸（*CONTACT_AUTOMATIC_NODES_TO_SURFACE），可以自動地識別粒子與表面的接觸，并防止粒子穿透殼單元，靜摩擦系數(shù)和動摩擦系數(shù)分別設置為0.5和0.3。

2.2 材料設置

環(huán)模材料是17CrNiMo，彈性模量E=2.1×1011，泊松比μ=0.25[14]。推板、料筒、?？椎炔考褂?0 號剛體材料模型（*MAT_RIGID），其中料筒和模孔部分為固定剛體，推板是可以單向移動的剛體。在運行過程中，推板以6 m/s的速度豎直向下壓制物料。

物料使用LS-DYNA中基于5號材料彈塑性本構模型（*MAT_SOIL_AND_FORM），結合已有研究中小麥粉的壓制曲線來定義模型中的力學參數(shù)，小麥粉應力應變曲線如圖5所示[15]，對應的壓制力位移曲線如圖6所示。

圖5 小麥粉應力應變曲線

圖6 小麥粉單軸壓制過程壓制力與位移曲線

以環(huán)模制粒機K15為例進行仿真，其相關參數(shù)為：壓輥半徑（r）為80 mm，環(huán)模半徑（R）為350 mm，壓輥與環(huán)模中心距（d）為92 mm，環(huán)模轉速（n）為300 r/min。根據(jù)公式（2）計算物料的壓縮速度，將其簡化為勻速運動，壓縮速度為6 m/s。因此，將推板下壓速度設置為6 m/s，每次下壓量為14 mm，留有1 mm 的空隙作為模輥間隙。

2.3 加載過程

在環(huán)模制粒機實際運行過程中，環(huán)模旋轉一周，壓輥僅將一小部分原材料壓入通道中。因此，物料進給是連續(xù)的[16]，所以在數(shù)值模擬時，每次壓制之前加入一定量的物料。

LS-DYNA的重啟動或者Dynain文件法可以記錄上一次計算的結果，并完全并入下一次計算中去。因而可以實現(xiàn)多次加料的過程，使用LS-PrePost 生成SPH粒子，第一次壓制時生成了4 992個SPH粒子，考慮到環(huán)模制粒過程中存在模輥間隙，之后每次壓制生成4 224個SPH粒子。第二加載物料如圖7所示。

圖7 加載物料

3 結果與分析

3.1 物料的流動情況

物料的流動情況可以由物料沿模孔軸向的流速直觀地體現(xiàn)出來，如圖8 所示。圖8 顯示了第一次壓制和第二壓制過程中物料流動速度的剖面圖，該圖顯示了物料由料筒進入?？字?，物料先是向中間流動，然后再進入模孔，這與Nielsen等[6-7]的試驗結果一致。Nielsen 等[7]使用紅色和藍色墨水分別對粉料進行著色，使擠壓后的顆粒形成不同顏色的分層：在非連續(xù)多次擠壓試驗中，每次擠壓前添加固定重量的染色粉料，紅色和藍色粉料交替加入。單孔擠壓試驗結果如圖9所示，通過比較數(shù)值模擬的擠壓過程流動情況和擠壓顆粒剖面圖，驗證了所提出的SPH模型。

圖8 第一次壓制和第二次壓制過程中的物料流動情況

圖9 單孔不連續(xù)擠壓試驗結果[7]

此外，倒角上方的穩(wěn)定原料的速度等值線在此區(qū)域中顯示出非常低的速度。表明在料筒底部形成楔形的積料區(qū)，積料區(qū)與?？椎菇沁B接到一起，且積料區(qū)的角度與?？椎菇堑慕嵌然疽恢隆Ｒ虼?，積料區(qū)是?？椎菇堑难由?，是物料在壓制過程中自動形成的，這與Nielsen等[6-7]的試驗結果一致。圖9中紅色圓圈處顯示了在實際單孔擠壓過程也會形成楔形積料區(qū)，且六次擠壓試驗結果顯示積料區(qū)的角度基本一致。六次擠壓中使用的物料相同，在倒角不同的模具中進行擠壓，表明楔形積料區(qū)的角度與倒角的角度沒有直接關系。

積料區(qū)的形成，有利于物料從料筒進入模孔，在工程應用中，環(huán)模?？椎牡菇切枰鶕?jù)物料的物理特性進行設計[17]，模擬試驗的結果顯示，小麥粉進行擠壓時宜選用?？椎菇菫?5°的環(huán)模。

3.2 物料的密度分布

在第一次和第二次擠壓結束后，物料的密度分布情況如圖10 所示。如圖10 顯示，物料密度沿?？纵S向從入口到出口不斷減小。第一次壓制后，物料密度的分布情況與流動情況相反，流速大的部分密度反而小。同時，與第一次壓制相比，第二次壓制后，物料的密度明顯增大。第二次壓制后，?？字虚g處的物料的密度大于邊緣處，與物料的流動情況較一致。

圖10 第一次壓制和第二次壓制結束時物料的密度分布情況

在單孔擠壓連續(xù)壓制過程中，需要背壓來啟動該過程。背壓是由材料在?？字蟹e聚產生的，這設定了壓力（預壓力）的要求，以克服?？變鹊哪Σ羀18]。第二次壓制與第一次壓制之間的差異，其原因可能是第一次壓制后?？兹肟谔幬锪系拿芏容^大，從而形成較大的背壓。

3.3 壓制次數(shù)與壓制力

為了探究壓制力與壓制次數(shù)的關系，進行了多次擠壓的單孔擠壓過程數(shù)值模擬，壓制過程中壓制力的變化情況如圖11 所示。從第二次壓制開始，壓制力與時間的關系曲線出現(xiàn)兩個峰，分析原因：上一次壓制完，料筒內物料的密度比模孔內物料的密度大，到達第1 個峰時，新加入的物料的密度達到前一次壓制后料筒內物料的密度；推板繼續(xù)下壓，新加入的物料與上次壓制后仍停留在料筒內的物料一起被壓縮，當推板的壓制力等于物料的黏結力及?？椎哪Σ亮r，出現(xiàn)第2 個峰，隨后，料筒內物料被擠入?？住?/p>

圖11 不同擠壓次數(shù)與壓制力

第5次壓制時，第1個峰并不明顯，且第2個峰峰值與第4 次基本一致，表明第4 次壓制完料筒內殘留的物料與?？變任锪系拿芏纫恢?，當新加入的物料與上次壓制后料筒內物料的密度一致后，物料的密度不再增加，而是直接被擠出。表明第4 次壓制后，進入穩(wěn)定擠壓狀態(tài)。

4 結論

提出了基于SPH的環(huán)模制粒成型數(shù)值模擬方法，以小麥粉為對象，結合K15型環(huán)模制粒機的結構參數(shù)進行了仿真分析。結果顯示，小麥粉在壓制過程中，物料會在料筒底部形成積料區(qū)，其作用與?？椎菇且粯涌梢允刮锪细玫亓鲃?；在進入穩(wěn)定擠壓狀態(tài)前，?？兹肟谔幬锪系拿芏雀?，形成更大的背壓，使下一次壓制時?？字虚g物料的密度更大；當進行4次壓制后，單孔擠壓進入穩(wěn)定擠壓狀態(tài)。仿真結果與已有研究結果一致，該方法適用于飼料粉體、秸稈顆粒等松散物料環(huán)模制粒過程的數(shù)值模擬，對其他領域松散物料模壓成型工藝的仿真有指導意義。