田志國(guó)， 王沫然

0 引 言

供油系統(tǒng)采用的多孔芯閥材料一般的制備過程是：顆粒物質(zhì)冷壓、燒結(jié).顆粒物質(zhì)的冷壓過程有很多方面的應(yīng)用，如藥片的制備、陶瓷工業(yè)、食品加工等.冷壓過程通常可以分為四步[1].第一步是由顆粒物質(zhì)的移動(dòng)和旋轉(zhuǎn)導(dǎo)致的重排過程，壓力較低，顆粒物質(zhì)尚沒有發(fā)生變形.第二步是隨著壓力升高，顆粒物質(zhì)之間互相鎖定，出現(xiàn)彈性變形，因此稱為彈性變形階段.隨著壓力的進(jìn)一步升高，顆粒物質(zhì)進(jìn)一步被擠壓而進(jìn)入塑性變形區(qū)域，這一步稱為塑性變形階段.最后，當(dāng)壓力足夠大時(shí)，顆粒物質(zhì)發(fā)生破裂.前人對(duì)這一過程進(jìn)行了很多實(shí)驗(yàn)研究，但實(shí)驗(yàn)研究通常比較耗時(shí)，并且實(shí)驗(yàn)研究通常是多因素混雜，難以研究單一因素的影響，因此，目前有很多研究者選擇采用數(shù)值模擬的方法對(duì)這一問題進(jìn)行研究，提出了一些連續(xù)介質(zhì)模型，如Walker模型[2]、Heckle模型[3]、Kawakita模型[4]等.這些連續(xù)介質(zhì)模型是宏觀模型，關(guān)注的是宏觀力學(xué)性質(zhì)與加載壓力之間的關(guān)系，忽略了內(nèi)部微結(jié)構(gòu)的演變過程.這類模型適用于研究冷壓結(jié)構(gòu)的宏觀性質(zhì)，但是在預(yù)測(cè)冷壓結(jié)構(gòu)的輸運(yùn)性質(zhì)時(shí)就不再適用.本研究主要關(guān)注冷壓結(jié)構(gòu)的流體輸運(yùn)性質(zhì)，其主要取決于內(nèi)部孔隙結(jié)構(gòu).因此，本研究需要采用微觀尺度或介觀尺度的數(shù)值模擬方法進(jìn)行研究.

離散單元法是一種介觀尺度方法，自1979年被Cundall[5]提出以來，已經(jīng)被廣泛用于顆粒系統(tǒng)的研究.離散單元法的基本思想是將材料視為離散顆粒的集合，通過計(jì)算施加在各個(gè)顆粒上的合力，并根據(jù)牛頓第二定律計(jì)算各個(gè)顆粒的速度和位移，進(jìn)而追蹤每個(gè)顆粒.離散單元法的核心是接觸力模型.本研究采用的接觸力模型是線性彈簧阻尼模型，是Cundall在1979年提出的.這個(gè)模型的優(yōu)點(diǎn)是足夠簡(jiǎn)單的同時(shí)又能定性地反映材料微觀結(jié)構(gòu)的變化.

研究流體輸運(yùn)性質(zhì)還需要進(jìn)行流動(dòng)的模擬.多孔介質(zhì)具有復(fù)雜的孔隙結(jié)構(gòu).在眾多的計(jì)算流體力學(xué)方法中，格子玻爾茲曼方法[6-9]具有易于并行計(jì)算和能夠高效處理復(fù)雜邊界條件等優(yōu)點(diǎn)，因此，本研究采用格子玻爾茲曼方法進(jìn)行流動(dòng)的模擬.

本研究結(jié)合離散單元法和格子玻爾茲曼方法來研究冷壓結(jié)構(gòu)的滲透率隨著壓力的變化.第二部分將介紹離散單元法和格子玻爾茲曼方法的基本概念，并介紹獲取冷壓結(jié)構(gòu)的方法和步驟.第三部分通過數(shù)值模擬得到等大球體面心立方堆積結(jié)構(gòu)的滲透率，并且模擬結(jié)果與理論計(jì)算的結(jié)果吻合，從而驗(yàn)證了格子玻爾茲曼方法的正確性.第四部分進(jìn)行參數(shù)敏感性分析，分析球體顆粒間摩擦角對(duì)冷壓結(jié)構(gòu)滲透率的影響.最后，給出研究結(jié)論.

1 數(shù)值模擬方法

1.1 離散單元法

以Cundall提出的線性彈簧阻尼模型為例.

兩個(gè)顆粒的法向接觸力的計(jì)算式為：

Fni=-Kniδnci-ζni(vci·nci)nci

(1)

切向接觸力的計(jì)算式為：

(2)

計(jì)算得到球體間的法向接觸力和切向接觸力后，利用牛頓第二定律計(jì)算各個(gè)顆粒的速度和角速度：

(3)

(4)

通過以上式子，在各個(gè)時(shí)間步上進(jìn)行顯式積分，可以完成離散單元法的求解過程.

1.2 格子玻爾茲曼方法

格子玻爾茲曼方法是一種基于流體密度分布函數(shù)對(duì)流動(dòng)進(jìn)行演化求解的數(shù)值模擬方法，其出發(fā)點(diǎn)是求解玻爾茲曼方程：

(5)

式中：c是粒子速度；f是分布函數(shù)，是粒子位置、速度和時(shí)間的函數(shù)；Ω(f)是碰撞項(xiàng)，通常是一個(gè)高度非線性的函數(shù).

對(duì)玻爾茲曼方程做Bhatnagar-Gross-Krook近似，并在格子空間進(jìn)行離散，得到單弛豫時(shí)間的格子玻爾茲曼方程為：

fi(x+ciΔt,t+Δt)=

(6)

(7)

式中：ωi是特定速度空間的權(quán)函數(shù)；c是格子速度；ρ是格子密度；u是格點(diǎn)上粒子的速度；ei為方向i上的單位向量.

可以證明，上式通過Chapman-Enskog展開可以回歸到Navier-Stokes方程[10-12].

1.3 使用離散單元法生成冷壓結(jié)構(gòu)

離散單元法發(fā)展至今，多款商業(yè)軟件和開源軟件已經(jīng)面世.本研究使用開源軟件，Yade進(jìn)行模擬[13]，使用的接觸力模型是線性彈簧阻尼模型，使用的接觸搜索算法是SAP，使用的積分算法是顯式積分.

生成冷壓結(jié)構(gòu)的步驟分為兩步：

第二步，是壓縮初始結(jié)構(gòu)到給定壓力.使用的邊界條件是伺服邊界條件，即監(jiān)測(cè)下壓面上的總壓力，當(dāng)總壓力大于給定壓力時(shí)，下壓面會(huì)適當(dāng)上移，以釋放壓力；當(dāng)總壓力小于給定壓力時(shí)，下壓面會(huì)適當(dāng)下移，以增大壓力，最終達(dá)到給定壓力值.

2 格子玻爾茲曼程序的驗(yàn)證

等大球體面心立方堆積結(jié)構(gòu)的滲透率存在解析解[14]：

(8)

式中：d是球體直徑，k是滲透率.

本研究的格子玻爾茲曼程序模擬得到等大球體面心立方堆積結(jié)構(gòu)的滲透率的值如圖1所示，橫軸N代表線分辨率，即每條邊有N個(gè)格點(diǎn)進(jìn)行解析.模擬結(jié)果顯示：在線分辨率較低時(shí)，滲透率的模擬值有較大的波動(dòng)；增大線分辨率可以使?jié)B透率的模擬值逐漸接近理論值.當(dāng)線分辨率N=50時(shí)，滲透率的模擬值與理論值的誤差小于5%，這意味著對(duì)于每個(gè)球體，用25×25×25的網(wǎng)格去解析，得到的模擬結(jié)果的誤差也會(huì)小于5%.在之后的模擬中，單個(gè)球體的解析度都會(huì)大于此值.

圖1 等大球體面心立方堆積結(jié)構(gòu)在不同分辨率下的模擬值與實(shí)驗(yàn)值Fig.1 Analytical results and simulation results of the face centered cubic (FCC) array under different resolution

3 參數(shù)敏感度分析

在球體顆粒冷壓結(jié)構(gòu)中，球體顆粒間的摩擦是影響冷壓結(jié)構(gòu)滲透率與壓力之間關(guān)系的重要參數(shù)之一.本研究通過改變球體顆粒間的摩擦角來研究摩擦對(duì)滲透率與壓力之間關(guān)系的影響.

在進(jìn)行模擬之前，首先需要確定代表體積單元(REV)的大小，在本研究中，即需要確定組成代表體積單元的最少顆粒數(shù)目.由文獻(xiàn)[15]中的研究可知，當(dāng)球體顆粒的數(shù)目達(dá)到4000時(shí)，生成過程的隨機(jī)性的影響就足夠小而可以忽略不計(jì)了.因此，本研究各個(gè)算例中的球體顆粒的數(shù)目均選擇為4000個(gè).

模擬中，離散單元法采用的參數(shù)為：楊氏模量：500 MPa；泊松比：0.1；顆粒密度：4000 kg/m3;重力加速度：9.81 m/s3，球體顆粒半徑分布滿足16～32 μm上的均勻分布，格子玻爾茲曼方法采用的參數(shù)為：網(wǎng)格尺寸：1 μm；流體密度：1000 kg/m3；運(yùn)動(dòng)粘性系數(shù)：10-6m2/s;邊界條件：無滑移邊界條件.

在冷壓結(jié)構(gòu)生成過程中，分別改變球體顆粒間的摩擦角為5°、10°、15°、30°，得到1 MPa、3 MPa、5 MPa、7 MPa、9 MPa、10 MPa、30 MPa、50 MPa、70 MPa、90 MPa對(duì)應(yīng)的冷壓結(jié)構(gòu)，再分別進(jìn)行滲流模擬，得到的模擬結(jié)果如圖2所示.

模擬結(jié)果顯示，在各個(gè)摩擦角下，冷壓結(jié)構(gòu)對(duì)應(yīng)的滲透率均隨著壓力的增大而減小，但相對(duì)滲透率的結(jié)果表明，不同的摩擦角下，冷壓結(jié)構(gòu)對(duì)壓力的響應(yīng)方式不一致：摩擦角越大，結(jié)構(gòu)的滲透率下降得越平緩，也就是結(jié)構(gòu)的穩(wěn)定性越好.文獻(xiàn)[15]研究了不同球體半徑分布的情況下，冷壓結(jié)構(gòu)滲透率隨壓力的變化，其結(jié)果顯示，在不同的球體半徑分布下，冷壓結(jié)構(gòu)的相對(duì)滲透率對(duì)壓力的響應(yīng)基本是一致的.這與本研究的結(jié)果不同，說明在冷壓過程中，球體顆粒間的摩擦角是一個(gè)需要著重考慮的參數(shù)，摩擦角部分決定了冷壓結(jié)構(gòu)滲透率對(duì)壓力的響應(yīng)方式.同時(shí)，本研究的結(jié)果還表明，在同樣的壓力下，摩擦角越大，冷壓結(jié)構(gòu)的滲透率越大，這與物理直觀判斷是一致的，進(jìn)一步表明本研究結(jié)果的合理性.

圖2 不同摩擦角下，冷壓結(jié)構(gòu)滲透率與壓力的關(guān)系Fig.2 The relationship between permeability and pressure under different friction angles

對(duì)模擬結(jié)果進(jìn)行擬合，采用文獻(xiàn)[15]中的擬合公式：

(9)

其中：k0是壓力為1 MPa對(duì)應(yīng)的滲透率，p0等于1 MPa，k是壓力p下對(duì)應(yīng)的滲透率.擬合公式中參數(shù)a和b的擬合結(jié)果如表1所示，結(jié)果顯示，在不同的摩擦角下，擬合公式的參數(shù)值會(huì)發(fā)生變化，參數(shù)a的相對(duì)變化程度較小，參數(shù)b則會(huì)出現(xiàn)較大的相對(duì)變化.

表1 擬合公式中參數(shù)a和b在不同摩擦角下對(duì)應(yīng)的值Tab.1 Values of parameters a and b in the fitting formula

4 結(jié) 論

本研究采用離散單元法和格子玻爾茲曼方法研究了球體顆粒堆積結(jié)構(gòu)受壓形變過程中滲透率的變化規(guī)律，著重討論了球體顆粒間摩擦角對(duì)結(jié)構(gòu)滲透率的影響.研究發(fā)現(xiàn)，摩擦角會(huì)改變冷壓結(jié)構(gòu)滲透率對(duì)壓力的響應(yīng)方式，在不同的摩擦角下，冷壓結(jié)構(gòu)相對(duì)滲透率與壓力的擬合曲線的參數(shù)值有較大的相對(duì)變化，并給出了在不同摩擦角下對(duì)應(yīng)的擬合公式.此研究結(jié)果對(duì)顆粒物質(zhì)冷壓制備過程有一定的指導(dǎo)意義：顆粒間的摩擦是一個(gè)需要考慮的參數(shù)，可以通過加入潤(rùn)滑劑等方式來調(diào)控制備成品的滲透率.比如聚酰亞胺顆粒，可以通過表面修飾來改變顆粒間摩擦角，以在更小的壓力下制備成品，節(jié)省資源.