馬洪斌 馬 巖 楊春梅

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

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

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

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

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

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

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

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

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

氣體在孔穴間的相互作用通?？梢圆捎肔enard-Jones相互作用勢來描述［6］，即

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

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

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

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

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

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

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

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

2 GCEMC模擬方法

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

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

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

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

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

采用VB編程產(chǎn)生隨機(jī)數(shù)，通過該程序產(chǎn)生任意多個隨機(jī)數(shù)，產(chǎn)生500個隨機(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é)模型為：

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

其取值的平均值為：

按照以下的步驟來搜索：

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

若 lk＜bk，則取：

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

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

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

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

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

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

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

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

GCEMC的分析過程如圖3所示。

圖3 GCEMC分析過程流程圖

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

熱磨時的溫度變化可以用等溫線描述，以柳杉為例，假設(shè)MDF的纖維單元是1 mm×1 mm×10 mm的理想形態(tài)，年產(chǎn)1.5萬m3MDF消耗原木2.7萬m3。2.7萬m3有纖維單元2.7萬億根。平均每秒鐘在磨盤內(nèi)有160 256根/m2纖維單元。而對于M200A型的熱磨機(jī)，年消耗木材14.4萬m3，平均每秒鐘在磨盤內(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Δ 時，細(xì)胞就裂解成纖維。此情況下在壓強(qiáng)很小的時候，相當(dāng)多的熱磨排氣體分子被孔穴吸入，而后隨著壓強(qiáng)的不斷增加，孔內(nèi)熱磨的氣體密度也會有一定的增加，但變化的輻度不大。此種孔徑的孔穴內(nèi)兩壁之間僅可以吸入一層熱磨氣體分子。熱磨氣體分子與孔壁之間的作用勢會非常小，這也就是說孔壁對熱磨排氣分子的吸引力會非常強(qiáng)，所以熱磨排氣分子會被很輕易地吸入到孔的內(nèi)部?？讖皆谶M(jìn)氣稍微增大一點(diǎn)后，排氣的過程馬上加快，細(xì)胞裂解在壓強(qiáng)很低時保持和接近為零。通過對每個細(xì)胞的分析就可以推算出整個機(jī)器的功率消耗。利用GCEMC仿真獲得的熱磨機(jī)功率消耗如表1所示?？芍捎貌煌吞柕臒崮C(jī)在不同的轉(zhuǎn)速下所消耗的功率。

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

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

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

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

4 結(jié)束語

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

有關(guān)人造板木片細(xì)胞熱磨模擬研究說明木片尺寸與木材細(xì)胞孔穴存在密切關(guān)系，在熱磨的情況下孔徑參數(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.