馬洪斌 馬 巖 楊春梅

(哈爾濱理工大學(xué)，哈爾濱，150080) (東北林業(yè)大學(xué))

熱磨法加工木纖維的工藝比較復(fù)雜，在造紙和中密度人造板(MDF)的纖維生產(chǎn)過(guò)程中包括了機(jī)械、物理和化學(xué)等作用，影響纖維質(zhì)量的因素較多，因此必須科學(xué)、合理地選擇工藝參數(shù)，才能生產(chǎn)出優(yōu)質(zhì)纖維［1-3］。目前，大型熱磨機(jī)的售價(jià)比較昂貴，這就限制了大型熱磨機(jī)實(shí)驗(yàn)的研究，使其無(wú)法像小型木工機(jī)械那樣，可以通過(guò)試驗(yàn)樣機(jī)來(lái)進(jìn)行樣機(jī)測(cè)試。因此，對(duì)大型熱磨機(jī)，經(jīng)常是通過(guò)計(jì)算機(jī)仿真來(lái)進(jìn)行模擬檢驗(yàn)，以得到最佳的工藝參數(shù)。

近代數(shù)學(xué)的最新理論成果是利用計(jì)算機(jī)模擬隨機(jī)現(xiàn)象產(chǎn)生隨機(jī)數(shù)據(jù)來(lái)進(jìn)行近似計(jì)算，這種方法叫做蒙特卡羅方法，也稱為隨機(jī)模擬(Random simulation)方法。巨正則蒙特卡羅方法(GCEMC，grand carlonical ensemble Monte Carlo method)是目前應(yīng)用比較普遍的方法，該法對(duì)于人造板熱磨裂解來(lái)說(shuō)是比較好的方法之一。筆者采用計(jì)算機(jī)模擬的GCEMC對(duì)MDF的熱磨工序進(jìn)行建模仿真，通過(guò)對(duì)人造板木片細(xì)胞熱磨的模擬研究，驗(yàn)證了在熱磨的情況下孔徑與排氣等溫線之間的關(guān)系。最后，通過(guò)柳杉熱磨工序加工實(shí)例，對(duì)筆者采用的方法進(jìn)行了驗(yàn)證。

1 人造板木片細(xì)胞孔穴模型

該方法是以分子間的相互作用即分子力為基礎(chǔ)，以統(tǒng)計(jì)方法或稱統(tǒng)計(jì)實(shí)驗(yàn)為手段，采用現(xiàn)代計(jì)算機(jī)對(duì)高溫高壓流體與木片間的高壓相互作用過(guò)程進(jìn)行模擬計(jì)算，從而給出在這種相互作用下流體系統(tǒng)達(dá)到平衡狀態(tài)時(shí)的有關(guān)性質(zhì)。筆者對(duì)GCEMC法用于人造板木片的微孔模型建模等做簡(jiǎn)要的介紹，并在此基礎(chǔ)上對(duì)人造板的多孔木片對(duì)氣體排氣作用的重要性做了研究。把高壓氣體作為分析木片結(jié)構(gòu)高壓裂解特性，即結(jié)構(gòu)表征的一種手段，對(duì)氣體分離與排氣引起的木材細(xì)胞孔內(nèi)的擴(kuò)膨脹造成人造板木片的孔穴劈裂作用進(jìn)行了分析。由于傳統(tǒng)理論對(duì)氣體分子在木片細(xì)胞微孔中的作用缺乏準(zhǔn)確的描述，因而在某些情況下得到的結(jié)果也就缺乏實(shí)際意義，特別是在微孔的情況下更是如此。在尋求對(duì)微孔排氣特性進(jìn)行準(zhǔn)確描述的過(guò)程中，最近出現(xiàn)的采用純統(tǒng)計(jì)力學(xué)的分子模擬以其獨(dú)特的優(yōu)勢(shì)受到人們的重視。

木片熱磨過(guò)程中氣體分子流動(dòng)模擬必須建立在有關(guān)人造板木片纖維之間的縫隙和木片細(xì)胞(通稱孔穴)孔穴的結(jié)構(gòu)模型，以及氣體分子與人造板孔穴之間的相互作用模型基礎(chǔ)上。隨著模擬研究工作的深入，有關(guān)孔結(jié)構(gòu)的模型也從簡(jiǎn)單到復(fù)雜。目前，人造板木片的孔穴分析并沒(méi)有真正的模型，為了分析方便，本研究采用的是應(yīng)用最多的圓形孔模型［4］，如圖1所示。木纖維細(xì)胞的六棱形的孔穴設(shè)為圓形孔，孔的分布的X、Y兩個(gè)方向均等，Z方向?yàn)橛邢揲L(zhǎng)的結(jié)構(gòu)。在X-Y平面h定義為孔直徑，各孔之間的距離近似相等。

圖1 不同孔徑下各孔穴相互作用勢(shì)的關(guān)系

由于實(shí)驗(yàn)所要進(jìn)行的分析是細(xì)胞級(jí)的研究，所以必須要區(qū)分細(xì)胞尺寸中直徑的不同定義，如考慮細(xì)胞孔壁的不同分布徑級(jí)和形狀使模型更符合實(shí)際。在今后的進(jìn)一步研究中，在相互作用模型中也可引入孔壁結(jié)構(gòu)的影響，以及孔穴體相互作用的影響等［5］，為了簡(jiǎn)化模擬計(jì)算，通常GCEMC模擬是以兩孔壁所在平面量得的孔直徑，即采用視頻檢測(cè)的方法確定細(xì)胞的孔徑。

從能量守恒原理可判定，單根木片的纖維束與單個(gè)氣體分子相互作用的勢(shì)能為零，所對(duì)應(yīng)的纖維束到氣體分子的距離對(duì)加壓蒸氣分子間相互作用的激撞直徑更接近于氣體分子所占據(jù)的真實(shí)直徑。在宏觀領(lǐng)域，這一差別可以略去不計(jì);而在孔穴的尺度范圍內(nèi)，則必須對(duì)這兩者加以區(qū)分，尤其是在進(jìn)行木片細(xì)胞孔徑分布分析時(shí)。在這樣一種圓形孔的細(xì)胞模型中，通常認(rèn)為孔壁表面與纖維束細(xì)胞表面是平行的，所以整個(gè)孔壁是纖維束表面沿Z方向平行排列而構(gòu)成。GCEMC是從分子間的相互作用出發(fā)討論流體與固體的相互作用。因此，作為GCEMC計(jì)算機(jī)模擬的一個(gè)輸入條件，必須已知?dú)怏w分子與孔壁之間相互作用的具體形式。在目前的計(jì)算機(jī)模擬中大多采用以Lenard-Jones作用勢(shì)為基礎(chǔ)的相互作用形式。

氣體在孔穴間的相互作用通常可以采用Lenard-Jones相互作用勢(shì)來(lái)描述［6］，即

式中：u表示兩個(gè)孔穴之間氣體熱膨脹的相互作用能;r表示兩孔穴之間的距離;εff為兩孔穴間相互熱量分子作用能的最小值，故也常稱之為相互作用勢(shì)阱的阱深;σff為兩孔穴之間氣體分子間相互作用勢(shì)為零所對(duì)應(yīng)的分子間距離，也稱碰撞直徑。σff是個(gè)變量，隨作用時(shí)間和氣體壓力變化而變化，開(kāi)始作用時(shí)σff為零，隨著作用的時(shí)間或壓力的加大，σff加大，當(dāng)孔穴在壓力下裂解σff和r相當(dāng)。

考慮到孔穴之間氣體的相互作用隨分子間的距離的增加而迅速衰減的具體情況，為了在計(jì)算機(jī)模擬的過(guò)程中節(jié)省機(jī)時(shí)，通常對(duì)式(1)進(jìn)行修正，得

式中：Rcut一般取5σff，一個(gè)細(xì)胞孔穴內(nèi)氣體分子與孔壁上另一個(gè)細(xì)胞孔穴內(nèi)分子間的相互作用仍適用關(guān)系式(1)，但有關(guān)參數(shù)做一些相應(yīng)的變化，得

式中，各量的物理意義與式(1)中相應(yīng)，只是這里表示兩個(gè)相鄰不同位置細(xì)胞的分子而已。如果采用圓形孔的模型，并將模擬具體化為人造板木片纖維細(xì)胞間的情況，則每個(gè)木片斜斷面可以看作是由纖維束無(wú)數(shù)細(xì)胞孔穴組成的無(wú)限大平面。

做進(jìn)一步的假設(shè)，認(rèn)為孔壁上的細(xì)胞孔穴是均勻分布的，不單是細(xì)胞個(gè)體之間的孔穴相互作用，相鄰孔穴之間又相互作用，不相鄰的r＜Rcut之內(nèi)的都對(duì)U(z)產(chǎn)生作用。細(xì)胞孔穴面密度(即單位面積上的細(xì)胞孔穴數(shù))為ρa(bǔ)，則一個(gè)孔壁和一個(gè)氣體分子的作用勢(shì)為

式中：k表示有相互作用氣體分子作用中心到被作用孔穴中心的距離;Δ為纖維束作用表面中心的間距。第1、第2項(xiàng)對(duì)應(yīng)孔壁中第一層纖維束層的貢獻(xiàn)，最后一項(xiàng)對(duì)應(yīng)第二細(xì)胞層及以后所有細(xì)胞層的貢獻(xiàn)。顯然式(4)隱含了一個(gè)假設(shè)，即孔壁是無(wú)限厚的。這與人造板木片細(xì)胞的真實(shí)分布情況是不符的。由于細(xì)胞管胞孔徑和木片的宏觀尺寸與細(xì)胞相比可以認(rèn)為是無(wú)窮大，作為一種對(duì)理想情況的研究是可以的。在圓形孔的假設(shè)下，每個(gè)孔有兩個(gè)孔距。因此，整個(gè)孔與孔的相互作用勢(shì)為

圖2所示細(xì)胞不同孔徑和相鄰細(xì)胞下各孔穴相互作用勢(shì)U與k的關(guān)系。可見(jiàn)，在不同的細(xì)胞孔徑下，總作用勢(shì)是式(2)限定以后的疊加。由于作用勢(shì)上的重疊程度不同，合成勢(shì)的曲線有明顯的差異。圖2所示孔徑越小，重迭越多，勢(shì)曲線相應(yīng)的阱深越大。特別是在很小的孔徑下，勢(shì)曲線的兩個(gè)極小合并為一個(gè)。一般的木片細(xì)胞的空穴不存在這種情況，這種勢(shì)曲線上的差異必定會(huì)導(dǎo)致不同孔徑下排氣特性有明顯的區(qū)別。

圖2 不同孔徑下各孔穴相互作用勢(shì)的關(guān)系

2 GCEMC模擬方法

微米級(jí)纖維板的纖維是將預(yù)熱處理后的木片在磨齒間進(jìn)行加工而形成的。在這里相對(duì)于木片本身來(lái)講，纖維在熱磨機(jī)內(nèi)所受到的作用力和作用次數(shù)一般是滿足下面兩個(gè)條件的樣本序列的馬爾科夫鏈：第一，處于磨盤(pán)上的木片均是一個(gè)有限的木片樣本集合，把這一集合稱之為木片的狀態(tài)空間;第二，在磨盤(pán)上除了第一個(gè)木片樣本點(diǎn)以外，其他任一個(gè)木片樣本點(diǎn)出現(xiàn)的幾率都僅與這個(gè)木片樣本前一次點(diǎn)所出現(xiàn)的位置相關(guān)。在馬爾科夫過(guò)程中，當(dāng)每一個(gè)木片達(dá)到過(guò)程順應(yīng)狀態(tài)，熱磨機(jī)中的木片在磨削過(guò)程中從其中的一個(gè)狀態(tài)開(kāi)始，經(jīng)過(guò)若干次后重新返回到這一狀態(tài)的概率為1時(shí)，那么這樣馬爾科夫鏈就會(huì)遵循于某個(gè)唯一的分布，即在磨削過(guò)程中，任一個(gè)初始狀態(tài)變化只要是遵循上述的馬爾科夫鏈，那么其中間不論經(jīng)歷的是怎樣的過(guò)程，最終都將逐步趨近于某個(gè)確定的狀態(tài)。如果滿足了馬爾科夫狀態(tài)，就可以實(shí)施蒙特卡羅模擬。

運(yùn)用蒙特卡羅方法模擬的宗旨是計(jì)算木片形成纖維過(guò)程切削力的數(shù)學(xué)期望值，根據(jù)大數(shù)法則和中心極限定理，采用隨機(jī)數(shù)序列對(duì)纖維的每一個(gè)狀態(tài)的抽樣值進(jìn)行狀態(tài)空間的建立。在這里，由于GCEMC是一個(gè)巨正則系統(tǒng)，所以在模擬過(guò)程中要保證分析系統(tǒng)的體積V、溫度T及系統(tǒng)的化學(xué)勢(shì)P不發(fā)生變化，這里的化學(xué)勢(shì)不變也等同于系統(tǒng)壓強(qiáng)P不變。模擬分析系統(tǒng)中所講的木片體積在這里指的就是模擬體的體積，一般取寬等于木片的寬度，X方向?yàn)槟酒?xì)胞端面最寬的方向，Y方向?yàn)槟酒?xì)胞端面最窄的方向同時(shí)也與X方向相垂直，木片可以看作是由一定長(zhǎng)度L的模擬胞所組成。隨著L的變大，將會(huì)增多處于模擬胞內(nèi)的氣體分子數(shù)，計(jì)算結(jié)果也會(huì)更加準(zhǔn)確，不足的是所需要的計(jì)算機(jī)時(shí)間會(huì)相應(yīng)地變長(zhǎng);反之，L變小后模擬胞內(nèi)的分子數(shù)就會(huì)變少，計(jì)算時(shí)間也會(huì)減少，這樣計(jì)算出來(lái)的統(tǒng)計(jì)誤差就會(huì)增加。在實(shí)際的模擬過(guò)程中，木片L的大小一般是要根據(jù)實(shí)際情況來(lái)選定。在研究木片裂解情況時(shí)，在給定的P、V、T下，要計(jì)算達(dá)到了裂解狀態(tài)時(shí)的模擬胞內(nèi)所含的氣體分子數(shù)，計(jì)算模擬相應(yīng)的細(xì)胞胞管裂解擴(kuò)張σff的變化，就可以得到裂解過(guò)程木片變成纖維的過(guò)程模擬。計(jì)算氣體膨脹的同時(shí)，還可以計(jì)算出與系統(tǒng)特性相關(guān)的一些其他的物理量，比如裂解木片的力學(xué)特性以及參數(shù)變化等等。

2.1 蒙特卡羅隨機(jī)數(shù)

用Monte Carlo方法優(yōu)化時(shí)，需要產(chǎn)生各種概率分布的隨機(jī)變量。對(duì)最簡(jiǎn)單、最基本的［0，1］上均勻分布隨機(jī)變量的抽樣值為隨機(jī)數(shù)。［0，1］上均勻分布的隨機(jī)變量的密度函數(shù)為：

數(shù)學(xué)期望為

采用VB編程產(chǎn)生隨機(jī)數(shù)，通過(guò)該程序產(chǎn)生任意多個(gè)隨機(jī)數(shù)，產(chǎn)生500個(gè)隨機(jī)數(shù)的分布的，隨機(jī)數(shù)均值為0.488 036 420371 87、方差為0.292 813 054 303 89、最大值和最小值分別為0.999 640 454 538 00和0.000 117 449 555 60。此方法產(chǎn)生的隨機(jī)數(shù)為偽隨機(jī)數(shù)，在實(shí)際應(yīng)用中，這些偽隨機(jī)數(shù)可以當(dāng)作真正的隨機(jī)數(shù)使用。

2.2 蒙特卡羅優(yōu)化仿真步驟

蒙特卡羅方法仿真，其數(shù)學(xué)模型為：

這是一個(gè)n維的不等式約束優(yōu)化問(wèn)題，設(shè)變量x=［x1，x2，…，xn］T上下界限為：

其取值的平均值為：

按照以下的步驟來(lái)搜索：

①在初值確定之后，要在設(shè)計(jì)變量的上下界限縮小的一個(gè)范圍內(nèi)開(kāi)始第j=1的第一輪隨機(jī)搜索i個(gè)點(diǎn)Xij。在這里，首先要確定縮小后的區(qū)間下限：

若 lk＜bk，則取：

再確定縮小后區(qū)間上限與下限之差：

若 lk+uk＞uk，則?。?/p>

其中，RND(x)為一隨機(jī)數(shù)，并且0＜RND(x)＜1。

②當(dāng)這個(gè)點(diǎn)得出之后，再用已知的約束條件來(lái)檢查其是否可行，若不可行，則需要重新產(chǎn)生一個(gè)新點(diǎn);如若可行，那么就計(jì)算其目標(biāo)函數(shù)f(x)，把計(jì)算出的目標(biāo)函數(shù)與事先給定的一個(gè)足夠大的數(shù)M相比較。若f(x)＞M，則不采用;若f(x)＜M，則可以采用，并且此時(shí)令f(x)?M。

③接下來(lái)檢查是否滿足迭代終止條件的要求，如若滿足就可以停止收縮，進(jìn)行步驟⑦;如若不然就需要進(jìn)行再一次輪的搜索。

⑤對(duì)每一次的搜索結(jié)果來(lái)講，不是新點(diǎn)不符合可行性或下降性的要求，就是得到了一個(gè)比上次更好的結(jié)果，如果搜索了i=imax次，就要進(jìn)一步圍繞此點(diǎn)縮小區(qū)間來(lái)進(jìn)行搜索，即令j+1?j。

⑥在當(dāng)搜索區(qū)間縮小jmax次之后，如果還沒(méi)有達(dá)到所需要求，就必須增加imax或jmax，重新進(jìn)行搜索。

⑦最后停止搜索，輸出結(jié)果。

GCEMC的分析過(guò)程如圖3所示。

圖3 GCEMC分析過(guò)程流程圖

3 GCEMC熱磨木片裂解的實(shí)例分析

熱磨時(shí)的溫度變化可以用等溫線描述，以柳杉為例，假設(shè)MDF的纖維單元是1 mm×1 mm×10 mm的理想形態(tài)，年產(chǎn)1.5萬(wàn)m3MDF消耗原木2.7萬(wàn)m3。2.7萬(wàn)m3有纖維單元2.7萬(wàn)億根。平均每秒鐘在磨盤(pán)內(nèi)有160 256根/m2纖維單元。而對(duì)于M200A型的熱磨機(jī)，年消耗木材14.4萬(wàn)m3，平均每秒鐘在磨盤(pán)內(nèi)有474 833根/m2纖維單元。它的細(xì)胞孔直徑為 Δ=0.037 mm，hmax=0.5Δ，hmin=0.15Δ，r=0.4Δ，σff=［0.15Δ，0.4Δ］，當(dāng) σff=0.4Δ 時(shí)，細(xì)胞就裂解成纖維。此情況下在壓強(qiáng)很小的時(shí)候，相當(dāng)多的熱磨排氣體分子被孔穴吸入，而后隨著壓強(qiáng)的不斷增加，孔內(nèi)熱磨的氣體密度也會(huì)有一定的增加，但變化的輻度不大。此種孔徑的孔穴內(nèi)兩壁之間僅可以吸入一層熱磨氣體分子。熱磨氣體分子與孔壁之間的作用勢(shì)會(huì)非常小，這也就是說(shuō)孔壁對(duì)熱磨排氣分子的吸引力會(huì)非常強(qiáng)，所以熱磨排氣分子會(huì)被很輕易地吸入到孔的內(nèi)部?？讖皆谶M(jìn)氣稍微增大一點(diǎn)后，排氣的過(guò)程馬上加快，細(xì)胞裂解在壓強(qiáng)很低時(shí)保持和接近為零。通過(guò)對(duì)每個(gè)細(xì)胞的分析就可以推算出整個(gè)機(jī)器的功率消耗。利用GCEMC仿真獲得的熱磨機(jī)功率消耗如表1所示。可知，采用不同型號(hào)的熱磨機(jī)在不同的轉(zhuǎn)速下所消耗的功率。

表1 熱磨機(jī)功率消耗

孔內(nèi)熱磨密度在壓強(qiáng)增大到了一定數(shù)值以后會(huì)迅速地增大。通常認(rèn)為，熱磨密度的迅速增大與孔內(nèi)兩壁表面的氣體壓力有關(guān)系。在壓強(qiáng)達(dá)到一定的數(shù)值之后，孔內(nèi)的熱磨密度會(huì)達(dá)到或接近于飽和狀態(tài)。同樣，隨著孔徑增大熱磨密度就會(huì)不斷增大，這說(shuō)明孔徑越大，孔的填充越大。有趣的是，當(dāng)存在于孔穴內(nèi)部的氣體達(dá)到裂解而破裂的時(shí)候，孔穴裂解所對(duì)應(yīng)著的壓強(qiáng)(在第二次熱磨排氣密度劇增時(shí)所對(duì)應(yīng)著的壓強(qiáng))也會(huì)隨著孔徑的不斷增大而增強(qiáng)，可是細(xì)胞的熱磨氣體分子密度相反則隨著孔徑的增大而不斷減小。這里除了最小的孔穴之外，其他熱磨氣體分子的飽和密度會(huì)隨著孔徑的增加而增大。細(xì)胞一旦破裂，氣體的密度和壓力馬上和環(huán)境平衡。從而形成纖維，纖維的尺寸如表2、表3所示。

表2 熱磨機(jī)尺寸變化

表3 熱磨機(jī)纖維數(shù)量變化

4 結(jié)束語(yǔ)

GCEMC計(jì)算機(jī)模擬方法，已成為一種研究人造板木片細(xì)胞特征的有力的數(shù)學(xué)和計(jì)算機(jī)模擬的工具。這個(gè)領(lǐng)域的研究現(xiàn)在已經(jīng)從開(kāi)始的理想化逐漸向?qū)嶋H材料結(jié)構(gòu)方面靠近［7］，已經(jīng)有相關(guān)研究來(lái)探討人造板孔壁壁厚及細(xì)胞結(jié)構(gòu)及其纖維分布的影響。模擬結(jié)果有的已用于分析實(shí)驗(yàn)得到的纖維生成過(guò)程的嘗試［8］?？偠灾?，GCEMC計(jì)算機(jī)模擬方法會(huì)在木片細(xì)胞研究方面發(fā)揮著更加積極的作用。

有關(guān)人造板木片細(xì)胞熱磨模擬研究說(shuō)明木片尺寸與木材細(xì)胞孔穴存在密切關(guān)系，在熱磨的情況下孔徑參數(shù)還會(huì)影響到排氣等溫線的一些基本特征，這對(duì)分析不同樹(shù)種的纖維得漿率具有重要意義。

［1］Widsten Petri.Oxidative activation of wood fibers for the manufacture of medium-density fiberborad［D］.Finland，Espoo：Helsinki University of Technology，2002.

［2］Sinou JJ，Thouverez F，Jezequel L.Non-linear analysis of a complex rotor/stator contact system［J］.Journal of Sound and Vibration，2004，278：1095-1129.

［3］Kim H，Cho M，Song S.Stability analysis of a turbine rotor system with Alford forces［J］.Journal of Sound and Vibration，2003，260(1)：167-182.

［4］Storace A F，Wisler D C，Shin H W，et al.Unsteady flow and whirl-inducing forces in axial-flow compressors：part I-experiment［J］.Journal of Turbomachinery，2001，123：433-445.

［5］馬巖.木材橫斷面六棱規(guī)則細(xì)胞數(shù)學(xué)描述理論研究［J］.生物數(shù)學(xué)學(xué)報(bào)，2002，17(1)：64-68.

［6］Jensen CR C，Papadopoulos G，Seaton NA，et al.Fundamentals of Adsorption［M］.Boston Massachusetts：Kluwer Academic Publishers，1996：422-439.

［7］Singh A P，Daniel G.The S2 layer in the Tracheid walls of Picea abies wood：Inhomogeneity in lignin Distribution and cell wall［J］.Holzforschung，2001，55(4)：373-378.

［8］Misato U W.Tangential Young’s Modulus of Coniferous Early Wood Investigated using Cell models［J］.Holzforschung，1999，53(2)：209-214.

［9］馬巖，楊春梅.蒙特卡羅方法在纖維熱磨中應(yīng)用的可行性初探［J］.林業(yè)科學(xué)，2006，42(10)：144-146.