馬洪斌，馬 巖，楊春梅

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

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

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

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

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

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

其取值的平均值為:

按以下步驟進(jìn)行搜索［4］:

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

若lk〈bk，則取:

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

若lk+uk〉uk，則取:

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

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

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

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

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

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

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

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

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

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

(1)Bs的確定。

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

(2)磨片每轉(zhuǎn)切斷長(zhǎng)的確定。

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

由下式可得磨片的每轉(zhuǎn)切斷長(zhǎng)為:

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

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

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

(5)刀齒齒長(zhǎng)的計(jì)算。

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

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

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

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

下面對(duì)切斷長(zhǎng)是否符合要求進(jìn)行驗(yàn)證。

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

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

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

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

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

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

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

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

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

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

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

齒高δ影響漿料在磨區(qū)的通過(guò)量和磨片的使用壽命，磨片齒角度β和熱磨機(jī)的泵送作用相關(guān)，磨片齒數(shù)Z和齒長(zhǎng)L是磨片的重要評(píng)價(jià)指標(biāo)，本研究主要選這四項(xiàng)參數(shù)進(jìn)行優(yōu)化。

(1)取設(shè)計(jì)變量為:

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

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

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

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

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

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

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

(4)確定約束條件。

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

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

上述就是動(dòng)盤、定盤的約束條件，公式 (13)和公式 (15)與已知纖維最小長(zhǎng)度有關(guān)，公式(14)與磨片的接觸面積有關(guān) (取一般經(jīng)驗(yàn)值)，引入到公式中并加以詳細(xì)的推導(dǎo)，在這里對(duì)動(dòng)、定盤的約束條件是一樣的，因?yàn)樗麄兪欠磳?duì)稱的。

3 蒙特卡羅隨機(jī)數(shù)產(chǎn)生

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

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

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

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

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

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

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

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

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

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

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

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

［2］李國(guó)云，王 忠.蒙特卡羅法在機(jī)械優(yōu)化設(shè)計(jì)中的應(yīng)用［J］.攀枝花學(xué)院學(xué)報(bào)，2009(6):63-66.

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

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

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

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

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

［8］史玉梅.圓環(huán)直齒熱磨機(jī)磨片的設(shè)計(jì)分析和結(jié)構(gòu)優(yōu)化［D］.長(zhǎng)春:吉林大學(xué)，2007.

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