纤维板热磨工序GCEMC法仿真1)

马洪斌 马 岩 杨春梅

(哈尔滨理工大学，哈尔滨，150080) (东北林业大学)

热磨法加工木纤维的工艺比较复杂，在造纸和中密度人造板(MDF)的纤维生产过程中包括了机械、物理和化学等作用，影响纤维质量的因素较多，因此必须科学、合理地选择工艺参数，才能生产出优质纤维［1-3］。目前，大型热磨机的售价比较昂贵，这就限制了大型热磨机实验的研究，使其无法像小型木工机械那样，可以通过试验样机来进行样机测试。因此，对大型热磨机，经常是通过计算机仿真来进行模拟检验，以得到最佳的工艺参数。

近代数学的最新理论成果是利用计算机模拟随机现象产生随机数据来进行近似计算，这种方法叫做蒙特卡罗方法，也称为随机模拟(Random simulation)方法。巨正则蒙特卡罗方法(GCEMC，grand carlonical ensemble Monte Carlo method)是目前应用比较普遍的方法，该法对于人造板热磨裂解来说是比较好的方法之一。笔者采用计算机模拟的GCEMC对MDF的热磨工序进行建模仿真，通过对人造板木片细胞热磨的模拟研究，验证了在热磨的情况下孔径与排气等温线之间的关系。最后，通过柳杉热磨工序加工实例，对笔者采用的方法进行了验证。

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

该方法是以分子间的相互作用即分子力为基础，以统计方法或称统计实验为手段，采用现代计算机对高温高压流体与木片间的高压相互作用过程进行模拟计算，从而给出在这种相互作用下流体系统达到平衡状态时的有关性质。笔者对GCEMC法用于人造板木片的微孔模型建模等做简要的介绍，并在此基础上对人造板的多孔木片对气体排气作用的重要性做了研究。把高压气体作为分析木片结构高压裂解特性，即结构表征的一种手段，对气体分离与排气引起的木材细胞孔内的扩膨胀造成人造板木片的孔穴劈裂作用进行了分析。由于传统理论对气体分子在木片细胞微孔中的作用缺乏准确的描述，因而在某些情况下得到的结果也就缺乏实际意义，特别是在微孔的情况下更是如此。在寻求对微孔排气特性进行准确描述的过程中，最近出现的采用纯统计力学的分子模拟以其独特的优势受到人们的重视。

木片热磨过程中气体分子流动模拟必须建立在有关人造板木片纤维之间的缝隙和木片细胞(通称孔穴)孔穴的结构模型，以及气体分子与人造板孔穴之间的相互作用模型基础上。随着模拟研究工作的深入，有关孔结构的模型也从简单到复杂。目前，人造板木片的孔穴分析并没有真正的模型，为了分析方便，本研究采用的是应用最多的圆形孔模型［4］，如图1所示。木纤维细胞的六棱形的孔穴设为圆形孔，孔的分布的X、Y两个方向均等，Z方向为有限长的结构。在X-Y平面h定义为孔直径，各孔之间的距离近似相等。

图1 不同孔径下各孔穴相互作用势的关系

由于实验所要进行的分析是细胞级的研究，所以必须要区分细胞尺寸中直径的不同定义，如考虑细胞孔壁的不同分布径级和形状使模型更符合实际。在今后的进一步研究中，在相互作用模型中也可引入孔壁结构的影响，以及孔穴体相互作用的影响等［5］，为了简化模拟计算，通常GCEMC模拟是以两孔壁所在平面量得的孔直径，即采用视频检测的方法确定细胞的孔径。

从能量守恒原理可判定，单根木片的纤维束与单个气体分子相互作用的势能为零，所对应的纤维束到气体分子的距离对加压蒸气分子间相互作用的激撞直径更接近于气体分子所占据的真实直径。在宏观领域，这一差别可以略去不计;而在孔穴的尺度范围内，则必须对这两者加以区分，尤其是在进行木片细胞孔径分布分析时。在这样一种圆形孔的细胞模型中，通常认为孔壁表面与纤维束细胞表面是平行的，所以整个孔壁是纤维束表面沿Z方向平行排列而构成。GCEMC是从分子间的相互作用出发讨论流体与固体的相互作用。因此，作为GCEMC计算机模拟的一个输入条件，必须已知气体分子与孔壁之间相互作用的具体形式。在目前的计算机模拟中大多采用以Lenard-Jones作用势为基础的相互作用形式。

气体在孔穴间的相互作用通常可以采用Lenard-Jones相互作用势来描述［6］，即

式中：u表示两个孔穴之间气体热膨胀的相互作用能;r表示两孔穴之间的距离;εff为两孔穴间相互热量分子作用能的最小值，故也常称之为相互作用势阱的阱深;σff为两孔穴之间气体分子间相互作用势为零所对应的分子间距离，也称碰撞直径。σff是个变量，随作用时间和气体压力变化而变化，开始作用时σff为零，随着作用的时间或压力的加大，σff加大，当孔穴在压力下裂解σff和r相当。

考虑到孔穴之间气体的相互作用随分子间的距离的增加而迅速衰减的具体情况，为了在计算机模拟的过程中节省机时，通常对式(1)进行修正，得

式中：Rcut一般取5σff，一个细胞孔穴内气体分子与孔壁上另一个细胞孔穴内分子间的相互作用仍适用关系式(1)，但有关参数做一些相应的变化，得

式中，各量的物理意义与式(1)中相应，只是这里表示两个相邻不同位置细胞的分子而已。如果采用圆形孔的模型，并将模拟具体化为人造板木片纤维细胞间的情况，则每个木片斜断面可以看作是由纤维束无数细胞孔穴组成的无限大平面。

做进一步的假设，认为孔壁上的细胞孔穴是均匀分布的，不单是细胞个体之间的孔穴相互作用，相邻孔穴之间又相互作用，不相邻的r＜Rcut之内的都对U(z)产生作用。细胞孔穴面密度(即单位面积上的细胞孔穴数)为ρa，则一个孔壁和一个气体分子的作用势为

式中：k表示有相互作用气体分子作用中心到被作用孔穴中心的距离;Δ为纤维束作用表面中心的间距。第1、第2项对应孔壁中第一层纤维束层的贡献，最后一项对应第二细胞层及以后所有细胞层的贡献。显然式(4)隐含了一个假设，即孔壁是无限厚的。这与人造板木片细胞的真实分布情况是不符的。由于细胞管胞孔径和木片的宏观尺寸与细胞相比可以认为是无穷大，作为一种对理想情况的研究是可以的。在圆形孔的假设下，每个孔有两个孔距。因此，整个孔与孔的相互作用势为

图2所示细胞不同孔径和相邻细胞下各孔穴相互作用势U与k的关系。可见，在不同的细胞孔径下，总作用势是式(2)限定以后的叠加。由于作用势上的重叠程度不同，合成势的曲线有明显的差异。图2所示孔径越小，重迭越多，势曲线相应的阱深越大。特别是在很小的孔径下，势曲线的两个极小合并为一个。一般的木片细胞的空穴不存在这种情况，这种势曲线上的差异必定会导致不同孔径下排气特性有明显的区别。

图2 不同孔径下各孔穴相互作用势的关系

2 GCEMC模拟方法

微米级纤维板的纤维是将预热处理后的木片在磨齿间进行加工而形成的。在这里相对于木片本身来讲，纤维在热磨机内所受到的作用力和作用次数一般是满足下面两个条件的样本序列的马尔科夫链：第一，处于磨盘上的木片均是一个有限的木片样本集合，把这一集合称之为木片的状态空间;第二，在磨盘上除了第一个木片样本点以外，其他任一个木片样本点出现的几率都仅与这个木片样本前一次点所出现的位置相关。在马尔科夫过程中，当每一个木片达到过程顺应状态，热磨机中的木片在磨削过程中从其中的一个状态开始，经过若干次后重新返回到这一状态的概率为1时，那么这样马尔科夫链就会遵循于某个唯一的分布，即在磨削过程中，任一个初始状态变化只要是遵循上述的马尔科夫链，那么其中间不论经历的是怎样的过程，最终都将逐步趋近于某个确定的状态。如果满足了马尔科夫状态，就可以实施蒙特卡罗模拟。

运用蒙特卡罗方法模拟的宗旨是计算木片形成纤维过程切削力的数学期望值，根据大数法则和中心极限定理，采用随机数序列对纤维的每一个状态的抽样值进行状态空间的建立。在这里，由于GCEMC是一个巨正则系统，所以在模拟过程中要保证分析系统的体积V、温度T及系统的化学势P不发生变化，这里的化学势不变也等同于系统压强P不变。模拟分析系统中所讲的木片体积在这里指的就是模拟体的体积，一般取宽等于木片的宽度，X方向为木片细胞端面最宽的方向，Y方向为木片细胞端面最窄的方向同时也与X方向相垂直，木片可以看作是由一定长度L的模拟胞所组成。随着L的变大，将会增多处于模拟胞内的气体分子数，计算结果也会更加准确，不足的是所需要的计算机时间会相应地变长;反之，L变小后模拟胞内的分子数就会变少，计算时间也会减少，这样计算出来的统计误差就会增加。在实际的模拟过程中，木片L的大小一般是要根据实际情况来选定。在研究木片裂解情况时，在给定的P、V、T下，要计算达到了裂解状态时的模拟胞内所含的气体分子数，计算模拟相应的细胞胞管裂解扩张σff的变化，就可以得到裂解过程木片变成纤维的过程模拟。计算气体膨胀的同时，还可以计算出与系统特性相关的一些其他的物理量，比如裂解木片的力学特性以及参数变化等等。

2.1 蒙特卡罗随机数

用Monte Carlo方法优化时，需要产生各种概率分布的随机变量。对最简单、最基本的［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 蒙特卡罗优化仿真步骤

蒙特卡罗方法仿真，其数学模型为：

这是一个n维的不等式约束优化问题，设变量x=［x1，x2，…，xn］T上下界限为：

其取值的平均值为：

按照以下的步骤来搜索：

①在初值确定之后，要在设计变量的上下界限缩小的一个范围内开始第j=1的第一轮随机搜索i个点Xij。在这里，首先要确定缩小后的区间下限：

若 lk＜bk，则取：

再确定缩小后区间上限与下限之差：

若 lk+uk＞uk，则取：

其中，RND(x)为一随机数，并且0＜RND(x)＜1。

②当这个点得出之后，再用已知的约束条件来检查其是否可行，若不可行，则需要重新产生一个新点;如若可行，那么就计算其目标函数f(x)，把计算出的目标函数与事先给定的一个足够大的数M相比较。若f(x)＞M，则不采用;若f(x)＜M，则可以采用，并且此时令f(x)⇒M。

③接下来检查是否满足迭代终止条件的要求，如若满足就可以停止收缩，进行步骤⑦;如若不然就需要进行再一次轮的搜索。

⑤对每一次的搜索结果来讲，不是新点不符合可行性或下降性的要求，就是得到了一个比上次更好的结果，如果搜索了i=imax次，就要进一步围绕此点缩小区间来进行搜索，即令j+1⇒j。

⑥在当搜索区间缩小jmax次之后，如果还没有达到所需要求，就必须增加imax或jmax，重新进行搜索。

⑦最后停止搜索，输出结果。

GCEMC的分析过程如图3所示。

图3 GCEMC分析过程流程图

3 GCEMC热磨木片裂解的实例分析

热磨时的温度变化可以用等温线描述，以柳杉为例，假设MDF的纤维单元是1 mm×1 mm×10 mm的理想形态，年产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Δ 时，细胞就裂解成纤维。此情况下在压强很小的时候，相当多的热磨排气体分子被孔穴吸入，而后随着压强的不断增加，孔内热磨的气体密度也会有一定的增加，但变化的辐度不大。此种孔径的孔穴内两壁之间仅可以吸入一层热磨气体分子。热磨气体分子与孔壁之间的作用势会非常小，这也就是说孔壁对热磨排气分子的吸引力会非常强，所以热磨排气分子会被很轻易地吸入到孔的内部。孔径在进气稍微增大一点后，排气的过程马上加快，细胞裂解在压强很低时保持和接近为零。通过对每个细胞的分析就可以推算出整个机器的功率消耗。利用GCEMC仿真获得的热磨机功率消耗如表1所示。可知，采用不同型号的热磨机在不同的转速下所消耗的功率。

表1 热磨机功率消耗

孔内热磨密度在压强增大到了一定数值以后会迅速地增大。通常认为，热磨密度的迅速增大与孔内两壁表面的气体压力有关系。在压强达到一定的数值之后，孔内的热磨密度会达到或接近于饱和状态。同样，随着孔径增大热磨密度就会不断增大，这说明孔径越大，孔的填充越大。有趣的是，当存在于孔穴内部的气体达到裂解而破裂的时候，孔穴裂解所对应着的压强(在第二次热磨排气密度剧增时所对应着的压强)也会随着孔径的不断增大而增强，可是细胞的热磨气体分子密度相反则随着孔径的增大而不断减小。这里除了最小的孔穴之外，其他热磨气体分子的饱和密度会随着孔径的增加而增大。细胞一旦破裂，气体的密度和压力马上和环境平衡。从而形成纤维，纤维的尺寸如表2、表3所示。

表2 热磨机尺寸变化

表3 热磨机纤维数量变化

4 结束语

GCEMC计算机模拟方法，已成为一种研究人造板木片细胞特征的有力的数学和计算机模拟的工具。这个领域的研究现在已经从开始的理想化逐渐向实际材料结构方面靠近［7］，已经有相关研究来探讨人造板孔壁壁厚及细胞结构及其纤维分布的影响。模拟结果有的已用于分析实验得到的纤维生成过程的尝试［8］。总而言之，GCEMC计算机模拟方法会在木片细胞研究方面发挥着更加积极的作用。

有关人造板木片细胞热磨模拟研究说明木片尺寸与木材细胞孔穴存在密切关系，在热磨的情况下孔径参数还会影响到排气等温线的一些基本特征，这对分析不同树种的纤维得浆率具有重要意义。

［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］马岩.木材横断面六棱规则细胞数学描述理论研究［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］.林业科学，2006，42(10)：144-146.