纖維板熱磨工序GCEMC法仿真1)

馬洪斌 馬 巖 楊春梅

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

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

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

1 人造板木片細胞孔穴模型

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

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

圖1 不同孔徑下各孔穴相互作用勢的關系

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

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

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

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

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

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

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

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

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

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

圖2 不同孔徑下各孔穴相互作用勢的關系

2 GCEMC模擬方法

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

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

2.1 蒙特卡羅隨機數

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

數學期望為

采用VB編程產生隨機數，通過該程序產生任意多個隨機數，產生500個隨機數的分布的，隨機數均值為0.488 036 420371 87、方差為0.292 813 054 303 89、最大值和最小值分別為0.999 640 454 538 00和0.000 117 449 555 60。此方法產生的隨機數為偽隨機數，在實際應用中，這些偽隨機數可以當作真正的隨機數使用。

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

蒙特卡羅方法仿真，其數學模型為：

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

其取值的平均值為：

按照以下的步驟來搜索：

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

若 lk＜bk，則?。?/p>

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

若 lk+uk＞uk，則取：

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

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

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

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

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

⑦最后停止搜索，輸出結果。

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

圖3 GCEMC分析過程流程圖

3 GCEMC熱磨木片裂解的實例分析

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

表1 熱磨機功率消耗

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

表2 熱磨機尺寸變化

表3 熱磨機纖維數量變化

4 結束語

GCEMC計算機模擬方法，已成為一種研究人造板木片細胞特征的有力的數學和計算機模擬的工具。這個領域的研究現(xiàn)在已經從開始的理想化逐漸向實際材料結構方面靠近［7］，已經有相關研究來探討人造板孔壁壁厚及細胞結構及其纖維分布的影響。模擬結果有的已用于分析實驗得到的纖維生成過程的嘗試［8］。總而言之，GCEMC計算機模擬方法會在木片細胞研究方面發(fā)揮著更加積極的作用。

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

［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ī)則細胞數學描述理論研究［J］.生物數學學報，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］馬巖，楊春梅.蒙特卡羅方法在纖維熱磨中應用的可行性初探［J］.林業(yè)科學，2006，42(10)：144-146.