馬洪斌，馬 巖，楊春梅

(1.哈爾濱理工大學機械動力工程學院，哈爾濱150080;2.東北林業(yè)大學研究生院，哈爾濱150040;3.東北林業(yè)大學林業(yè)與木工機械工程技術中心，哈爾濱150040)

蒙特卡羅方法 (Monte Carlo Method)被廣泛應用在數(shù)學、物理、工程技術以及管理等諸多方面［1］。這種方法的特點是如果模擬的次數(shù)越多，那么其計算出的結果可靠性就會越大。蒙特卡羅方法特別適用于在計算機上對某些大型項目和新產品進行開發(fā)。針對一些含有大量不確定因素的復雜決策系統(tǒng)進行風險模擬分析和問題優(yōu)化蒙特卡羅方法也同樣適用。采取這種方法還可用于機械優(yōu)化設計中有關隨機方向搜索問題的研究［2，3］。

關于隨機方向搜索優(yōu)化的方法在各種參考文獻中已有詳細的介紹［4］。蒙特卡羅方法用于隨機方向搜索方面有以下的特點:可以作為無約束優(yōu)化方法，也可作為約束優(yōu)化方法，再者其程序短、效率高、概念清晰。文中研究了蒙特卡羅方法的原理和其算法框圖，并用其對熱磨機磨片齒形進行結構優(yōu)化設計。

1 蒙特卡羅優(yōu)化算法原理

若有一優(yōu)化問題，其數(shù)學模型為［5-7］:

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

其取值的平均值為:

按以下步驟進行搜索［4］:

(1)在初值給定以后，我們依據(jù)設計變量所確定的上、下界限來縮小一個范圍，開始進行第j=1的第一輪隨機搜索，得到i個點Xij。第一個確定了的縮小后的區(qū)間下限為:

若lk〈bk，則取:

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

若lk+uk〉uk，則取:

縮小后的區(qū)間定后，就可求到一個點Xijk(i=1～500，j=1～8，k=1～n):

RND(x)為一隨機數(shù)，并且0〈RND(x)〈1。

(2)當獲得此點后，再用約束條件來檢查其是否可行，若不可行則再產生一個新點;若可行就以此計算出其目標函數(shù)f(x)，把這個目標函數(shù)與事先給定的一個足夠大的數(shù)M相比較，此時，若f(x)〉M則不采用，若f(x)〈M則采用，并且令f(x)?M。

(3)查看是否滿足迭代終止條件的要求，若滿足就可以停止收縮，然后轉到步驟⑦;若不滿足則再一次進行新的一輪搜索。

(5)針對每次搜索的結果，要不就是新點不符合可行性或下降性的要求，再有就是會都得到比上一次更好的結果，在搜索i=imax次以后，我們可以進一步圍繞此點來縮小區(qū)間開始新的搜索，令j+1?j。

(6)在當搜索區(qū)間被縮小jmax次之后，如果仍沒達到要求，我們可以把imax或jmax增加，然后重新搜索。

(7)搜索停止，然后輸出結果。

2 蒙特卡羅熱磨機磨片優(yōu)化數(shù)學模型

本文所要設計磨片所要處理的漿種是漂白硫酸鹽針葉材，這種原料的纖維平均長度為2.3～3.8 mm，平均寬度為39.8～47μm。采用的打漿方式為中等纖維長度的半粘狀、半游離狀漿。所配用電機功率為1400 kW，其轉速為1500 r/min。磨片的規(guī)格為:外徑1 070 mm，內徑400 mm。選用材料為00Cr27Mo，該種材料的各項性能參數(shù)為:屈服極限為645 MPa，強度極限為810 MPa，許用應力為712 MPa，楊氏模量為805 GPa，泊松比為0.32，密度為 8 315 kg/m3。

2.1 未優(yōu)化磨片參數(shù)的選定

(1)Bs的確定。

處理的漿種為半粘狀、半游離狀漂白硫酸鹽針葉木漿，齒刃比負荷 Bs取2.5～3.0 J/m之間的2.6 J/m。

(2)磨片每轉切斷長的確定。

齒刃比負荷 Bs已選定。配用電機功率為1 400 kW，轉速為1 500 r/min。該機的空載功率一般為240 kW，考慮實際生產中用于纖維分離的總功率一般在70%額定功率和滿載之間，因此取中間值為1190 kW，則用于分離纖維的有效功率為:

由下式可得磨片的每轉切斷長為:

(3)磨片齒紋類型的選擇。

本文研究的磨片是圓環(huán)分區(qū)直長齒類型的磨片。該磨盤圓環(huán)被劃分為破碎區(qū)、粗磨區(qū)和精磨區(qū)3個圓環(huán)區(qū)。各環(huán)的長度取值為:破碎區(qū)取90 mm，粗磨區(qū)取70 mm，精磨區(qū)取175 mm。粗磨區(qū)的刀齒由精磨區(qū)的刀齒相間延長，而破碎區(qū)的刀齒則由粗磨區(qū)的刀齒相間延長直至內圓邊，即粗磨區(qū)的刀齒數(shù)為精磨區(qū)的刀齒數(shù)的一半，而破碎區(qū)的刀齒數(shù)為粗磨區(qū)的刀齒數(shù)的一半，這樣取值便于設計和加工［8］。

(4)刀齒角度的選擇。本文中刀齒角度取10°。

(5)刀齒齒長的計算。

首先計算精磨區(qū)齒長，磨片半徑 (精磨區(qū)外徑)為535 mm，粗磨區(qū)外徑535-175=360 mm，刀齒角度10°。因此，由余弦定理可以算得:精磨區(qū)齒長為179.06 mm。以同樣的方法可以算得粗磨區(qū)和破碎區(qū)的長度分別為71.53 mm和92.04 mm。

(6)齒寬、槽寬和齒數(shù)的確定。

精磨區(qū)的齒寬取4 mm，槽寬取8.5 mm，那么一個齒節(jié)長就是12.5 mm。由于磨片的齒數(shù)一般較多，所以齒節(jié)長與齒周節(jié)長可以被近似看作相等為12.5 mm。那么精磨區(qū)的齒數(shù)就為1 070×3.14/12.5=268.78，在這里為了便于等分，取齒數(shù)為268。

由于精磨區(qū)的齒數(shù)為268個，故對應的粗磨區(qū)的齒數(shù)為268/2=134個，破碎區(qū)的齒數(shù)為134/2=67個。那么，粗磨區(qū)外圓邊一個齒節(jié)長約為(1 070-350) × 3.14/134=16.8 mm，若粗磨區(qū)的齒寬仍取4 mm，則粗磨區(qū)的槽寬就為12.8 mm。同理，破碎區(qū)上單個齒節(jié)長為 (1 070-490)×3.14/67=27.2 mm，如果破碎區(qū)的齒寬取8 mm，則槽寬為19.2 mm。

下面對切斷長是否符合要求進行驗證。

已知精磨區(qū)齒長為179.06 mm，齒數(shù)為268個，則精磨區(qū)的切斷長為:

268 ×268 ×0.17 906=12 860.80 544 m/r。

粗磨區(qū)齒長為71.53 mm，齒數(shù)為134個，則粗磨區(qū)的切斷長為:

134 × l34 ×0.07 153=1 284.3 927 m/r。

破碎區(qū)齒長為92.04 mm，齒數(shù)為67個，則破碎區(qū)的切斷長為:

67 ×67 ×0.09 204=413.1 677 m/r。

所以，總的切斷長為:14 558.36 568 m/r。該值與需要切斷長的偏差為0.39%，基本上吻合要求。因此，各環(huán)的齒寬、槽寬、齒數(shù)確定依次為:精磨區(qū)的齒寬為4 mm，槽寬為8.5 mm，齒數(shù)為268個;粗磨區(qū)的齒寬為4 mm，槽寬為12.8 mm，齒數(shù)為134個;破碎區(qū)的齒寬為8 mm，槽寬為19.2 mm，齒數(shù)為67個。

(7)齒高的確定。本文中齒高取7 mm。

(8)數(shù)據(jù)匯總。齒紋類型及刀齒數(shù)據(jù)匯總見表1所示。

表1 齒紋類型及刀齒數(shù)據(jù)Tab.1 Types of tooth line and data of disc tooth

2.2 磨片參數(shù)優(yōu)化模型

齒高δ影響漿料在磨區(qū)的通過量和磨片的使用壽命，磨片齒角度β和熱磨機的泵送作用相關，磨片齒數(shù)Z和齒長L是磨片的重要評價指標，本研究主要選這四項參數(shù)進行優(yōu)化。

(1)取設計變量為:

公式 (8)中表示磨盤內動盤和定盤刀齒的數(shù)目，其是參與熱磨的最重要參數(shù);L表示磨盤第i段內每秒切斷長，它往往決定了纖維的長度;δ分別是動盤刀齒的高度或定盤刀齒的高度，具體數(shù)值由磨盤的強度來確定;β是動盤和定盤上刀齒和半徑方向所成的夾角。整個數(shù)學模型的建立就是要對所設的參數(shù)和設計變量間的對應關系進行蒙特卡羅的優(yōu)化，以便后續(xù)推理論證的實現(xiàn)。

(2)其目標函數(shù)為:

由公式 (10)、公式 (11)和公式 (12)以及破碎區(qū)的取值90 mm，粗磨區(qū)的取值70 mm，精磨區(qū)的取值175 mm，再者粗磨區(qū)的刀齒數(shù)為精磨區(qū)的刀齒數(shù)一半，破碎區(qū)的刀齒數(shù)為粗磨區(qū)的刀齒數(shù)一半，可以得到上述公式 (9)所列的目標函數(shù)。由于V=πDn及Np均已知，磨片動盤的破碎區(qū)、粗磨區(qū)、精磨區(qū)每個區(qū)內齒數(shù)Z、齒高L、齒長δ和定盤的相同，β大小相等方向相反，式中引入的這些設計變量做為優(yōu)化的目標函數(shù)，這些設計變量與所需論證的目標函數(shù)的關系可以保證目標函數(shù)的齒刃比負荷為最小值。Z為磨盤動盤和定盤刀齒數(shù)目;L為刀齒長度，m;δ為分別為動盤刀齒的高度或定盤刀齒的高度，m;β為動盤和定盤上的刀齒和半徑方向所成的夾角;具體的上下界限數(shù)值在約束方程中給出，經過蒙特卡羅的推理與本文的基本假設，作為蒙特卡羅數(shù)學模型推導出來的優(yōu)化目標函數(shù)是實現(xiàn)最佳求解的前提。

式中:Fi為磨盤第i段內刀齒的接觸面積，m2;Zpi，Zci為磨盤第i段內動盤刀齒數(shù)和定盤刀齒數(shù);δpi，δci為磨盤第i段內動盤刀齒和定盤刀齒的厚度，m;Li為磨盤第i段刀齒的長度，m;βpi，βci為動盤和定盤上刀齒和圓盤直徑方向所成角度。

式中:Np為纖維分離時的有用功率，W;μ為纖維分離系數(shù);V為動盤圓周速度，m/s;F為動盤刀齒和定盤刀齒的接觸面積，m2。

式中:Pa為a區(qū)內刀齒間的平均比壓，Pa。

(3)確定設計變量的上下界限作為輸入?yún)?shù):

(4)確定約束條件。

動盤和定盤第i區(qū)約束條件為:

式中:Z為磨盤第i段內動盤和定盤刀齒數(shù)目;n為熱磨機動盤轉速，r/sec;L為第i段內動盤刀齒長度和定盤刀齒長度，m;δ為分別為第i段內動盤刀齒高度和定盤刀齒高度，m;β為動盤和定盤上的刀齒和半徑方向所成的角度;

上述就是動盤、定盤的約束條件，公式 (13)和公式 (15)與已知纖維最小長度有關，公式(14)與磨片的接觸面積有關 (取一般經驗值)，引入到公式中并加以詳細的推導，在這里對動、定盤的約束條件是一樣的，因為他們是反對稱的。

3 蒙特卡羅隨機數(shù)產生

在用蒙特卡羅方法進行優(yōu)化的時候，要產生各種概率分布的隨機變量。本研究取基本的 ［0，1］上均勻分布的隨機變量的抽樣值作為隨機數(shù)［9］。

為了產生隨機數(shù)，本文在EXCEL中采用VB編程產生隨機數(shù)，產生500個隨機數(shù)的分布如圖1所示，隨機數(shù)均值為0.48803642037187、方差為0.29281305430389、最大值和最小值分別為0.99964045453800 和 0.00011744955560。此方法產生的隨機數(shù)為偽隨機數(shù)，在實際應用中這些偽隨機數(shù)可以當作真正的隨機數(shù)使用。如圖1所示。

圖1 500個隨機數(shù)分布圖Fig.1 Distribution of 500 random number

4 蒙特卡羅優(yōu)化算法框圖

根據(jù)上述的蒙特卡羅算法原理，可以總結畫出蒙特卡羅優(yōu)化算法的框圖，如圖2所示。

圖2 蒙特卡羅優(yōu)化算法框圖Fig.2 Chart of optimization algorithm for Monte Carlo

用Visual Basic編制蒙特卡羅算法程序，可以得到經過優(yōu)化后的精磨區(qū)齒角度、齒長、齒數(shù)、齒高的數(shù)據(jù)，通過進一步計算可以得到粗磨區(qū)和破碎區(qū)的相應數(shù)據(jù)，具體數(shù)據(jù)見表2，定盤的參數(shù)與動盤的參數(shù)值相同，齒角度取相反值。

表2 優(yōu)化后的齒紋類型及刀齒數(shù)據(jù)Tab.2 Types of tooth line and data of disc tooth after optimization

通過表1和表2的數(shù)據(jù)可知，該設計外徑42英寸的磨片在經過優(yōu)化后的齒角度增加了3°，精磨區(qū)、粗磨區(qū)以及破碎區(qū)的齒長和齒高都有所增加，而齒數(shù)則有所減少，在此條件下對纖維的處理時間會比較理想，但如果再增加磨片尺寸，那么相應的每秒切斷長和平均比壓都會隨之升高，熱磨機對纖維的處理時間也就會隨之增長。表2中提到的“優(yōu)化后的齒紋類型及刀齒數(shù)據(jù)”在上述的論證中并沒有被特定的指出，但在實際工程中，齒紋類型是圓環(huán)分區(qū)的形狀，本文采用的是最簡單的直長齒。刀齒數(shù)據(jù)包括表中的所列數(shù)據(jù)。鑒于對優(yōu)化模型的精度要求，使得文中的理論數(shù)據(jù)在應用中要略有修改。

5 結束語

本文首先介紹了蒙特卡羅方法的算法和原理，并針對傳統(tǒng)設計方法的局限性，采用蒙特卡羅方法來進行熱磨機磨片磨齒部分的優(yōu)化設計。建立了蒙特卡羅熱磨機磨片優(yōu)化設計數(shù)學模型，并基于此模型對熱磨機磨片進行了優(yōu)化設計。在熱磨機磨片優(yōu)化設計過程中引入了蒙特卡羅方法。依據(jù)確定出的熱磨機磨片特征參數(shù)，通過產生隨機數(shù)并建立蒙特卡羅數(shù)學模型求解，從而獲得優(yōu)化了的磨片參數(shù)值。

［1］黃宏新，曹澤星.優(yōu)化量子Monte Carlo計算中的波函數(shù)［J］.湖南大學師范大學自然科學學報，1996(4):62-65.

［2］李國云，王 忠.蒙特卡羅法在機械優(yōu)化設計中的應用［J］.攀枝花學院學報，2009(6):63-66.

［3］丁柏群，張 鵬.蒙特卡羅模擬信號交叉口的服務水平［J］.森林工程，2006(3):35-37.

［4］劉星榮，方 珉.蒙特卡羅法在汽車變速器可靠性優(yōu)化設計中的應用［J］.輕型汽車技術，1998(5):4-8.

［5］覃 嶺，周 華，謝君生，等.汽車自動變速器檢測設備傳動系統(tǒng)的優(yōu)化設計［J］.現(xiàn)代機械，2008(3):9-12.

［6］張利偉，王 偉.基于參數(shù)化建模有限元優(yōu)化技術［J］.水科學與工程技術，2008(2):18-20.

［7］李宜海.機械優(yōu)化設計的數(shù)學建模及求解中的幾個問題［J］.煤礦現(xiàn)代化，2008(3):80-81.

［8］史玉梅.圓環(huán)直齒熱磨機磨片的設計分析和結構優(yōu)化［D］.長春:吉林大學，2007.

［9］姜 慧，崔國民，倪 錦.能量系統(tǒng)廣義換熱網絡優(yōu)化的蒙特卡羅SCDD法［J］.華北電力大學學報:自然科學版，2010(1):92-95.