蒙特卡羅法在熱磨機磨片優化設計中的應用

馬洪斌，馬 巖，楊春梅

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

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

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

1 蒙特卡羅優化算法原理

若有一優化問題，其數學模型為［5-7］:

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

其取值的平均值為:

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

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

若lk〈bk，則取:

再確定縮小后區間上限與下限之差:

若lk+uk〉uk，則取:

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

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

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

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

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

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

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

2 蒙特卡羅熱磨機磨片優化數學模型

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

2.1 未優化磨片參數的選定

(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)磨片齒紋類型的選擇。

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

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

(5)刀齒齒長的計算。

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

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

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

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

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

已知精磨區齒長為179.06 mm，齒數為268個，則精磨區的切斷長為:

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

粗磨區齒長為71.53 mm，齒數為134個，則粗磨區的切斷長為:

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

破碎區齒長為92.04 mm，齒數為67個，則破碎區的切斷長為:

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

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

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

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

表1 齒紋類型及刀齒數據Tab.1 Types of tooth line and data of disc tooth

2.2 磨片參數優化模型

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

(1)取設計變量為:

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

(2)其目標函數為:

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

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

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

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

(3)確定設計變量的上下界限作為輸入參數:

(4)確定約束條件。

動盤和定盤第i區約束條件為:

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

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

3 蒙特卡羅隨機數產生

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

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

圖1 500個隨機數分布圖Fig.1 Distribution of 500 random number

4 蒙特卡羅優化算法框圖

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

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

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

表2 優化后的齒紋類型及刀齒數據Tab.2 Types of tooth line and data of disc tooth after optimization

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

5 結束語

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

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

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

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

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

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

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

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

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

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