基于蒙特卡洛模擬與響應面方法的公差建模

余治民　劉子建　董思科　李斯明　艾彥迪

湖南大學汽車車身先進設計制造國家重點實驗室,長沙,410082

基于蒙特卡洛模擬與響應面方法的公差建模

余治民劉子建董思科李斯明艾彥迪

湖南大學汽車車身先進設計制造國家重點實驗室,長沙,410082

針對現有機床精度設計方法實用性不強的問題,提出了一種基于蒙特卡洛模擬與響應面方法的公差建模方法。采用基于數學定義的公差分析理論建立公差的變動不等式與約束不等式;運用蒙特卡洛模擬法進行仿真試驗，模擬實際公差表面的變動，生成公差變動要素的實際變動區間;以公差與試驗得到的公差變動要素實際變動區間帶寬值為建模樣本，運用響應面方法建立兩者間的響應面模型;對一個典型實例進行分析,分析結果表明,該方法符合工程實際，具有較高的建模精度及技術經濟性。

蒙特卡洛模擬；響應面法；小位移旋量；約束不等式；變動不等式

0　引言

科學技術的進步促使機床朝著高速和高精度方向發展,機床精度性能的設計與提高變得日益重要。國內外學者在機床精度設計領域開展了大量的研究,取得了一定的進展。文獻[1-3]運用多體運動學理論建立了機床的誤差傳遞模型;文獻[4-5]分析了加載時工作零件的變形量,將其轉變為雅可比旋量修正量, 通過對雅可比旋量公差模型在實際工況下的修正, 建立了基于雅可比旋量和實際工況的裝配體公差數學模型;文獻[6]給出了三維公差累積運動學模型的一般表示式,研究了該模型在公差優化設計中的應用。但上述研究都將機床零部件間的幾何變量約束統一用6項公差變動要素(3項平動誤差與3項轉動誤差)進行描述,而公差變動要素與零部件公差之間的關系并不明確,無法建立反映零部件公差與機床精度映射關系的數學模型,也就無法有效地對機床零部件公差進行合理分配，只能根據并不完善的公差模型來進行公差分配,可靠性低。因此,建立能用于工程實際且定義嚴格的公差模型，對機床的精度分析與分配研究具有重要意義。

公差建模的關鍵是可以對滿足公差的公差變動要素作出正確的解釋[7],即公差變動要素如何在公差域中變動。文獻[8-10]運用基于數學定義的公差建模方法，用變動不等式和約束不等式嚴謹地描述了公差變動要素與公差間的關系，建立了不同類型公差的數學模型,但由于約束不等式中變動要素存在變動順序不確定性,難以得到公差變動要素實際變動區間與公差間的具體函數關系,因此該模型無法用于直接指導精度設計。

本文以平面尺寸公差與平面度公差為研究對象，在深入研究基于數學定義的公差分析方法的基礎上,將蒙特卡洛模擬與響應面方法應用于公差建模，考慮公差原則及約束條件等因素對公差變動要素的影響，建立變動要素實際變動區間帶寬與公差間的響應面模型，以解決現有精度設計方法實用性不強的問題，并通過典型實例驗證該方法的有效性。

1　基于小位移旋量的公差數學表示

三維空間中，物體表面可抽象為點、線、面等基本要素,并有3個沿坐標軸方向的平動自由度和3個繞坐標軸的轉動自由度。用d=(u,v,w)和θ=(α,β,δ)分別表示平動與轉動自由度的微小變動矢量，則這兩組矢量組成的合成矢量D=(θ,d)=(α,β,δ,u,v,w)稱為小位移旋量(smalldisplacementtorsors,SDT),α、β、δ、u、v、w為SDT的旋量參數。

公差域表示公差實際表面脫離名義表面變動的范圍或區域，而物體表面的點、線、面等基本要素的變動量均可以通過各自的SDT表示。據此，Bourdet 在1996年首次將SDT引入到公差領域，提出了基于SDT的公差數學表示方法[11],即用旋量參數的變動來描述公差實際表面在公差域中的變動。旋量參數也稱為公差變動要素。表1列出了幾種形位公差對應的公差域和SDT。

2　基于蒙特卡洛模擬與響應面方法的公差建模

公差分為尺寸公差、形狀公差和位置公差。針對尺寸公差和位置公差，分析時可以假設變動后公差實際表面的形狀保持不變，即平面變動后依舊為平面，直線變動后依舊為直線，但對于形狀公差，這一假設并不成立[12]。為此，本文從兩類公差中各取一例，詳細闡述基于蒙特卡洛模擬與響應面方法的公差建模步驟,其他類型的公差建?？梢来祟愅?。

2.1平面尺寸公差建模

2.1.1求解公差變動不等式與約束不等式

平面尺寸公差規定的公差域形式如圖1所示。圖中,TU、TL分別為上下偏差,公差T=TU+TL;D為基本尺寸;局部坐標系的z軸與平面法向平行。對于矩形平面,坐標系原點在平面中心位置;對于不規則平面,可用一等效矩形替代,則原點在等效矩形的中心位置。z=0表示名義公差平面;z(x,y)表示實際變動平面。由表1可知,變動平面z(x,y)對應的SDT為(α,β,0,0,0,w)T。尺寸公差用數學表示如下:

-TL≤z(x,y)≤TU

(1)

-(TL+TU)≤Δz(x,y)≤TL+TU

(2)

表1　不同類型公差及其對應的公差域和SDT

圖1　平面尺寸公差域

其中,z(x,y)=dz+xβ+yα為實際變動平面方程,dz為實際變動平面中心點與參考坐標系原點在z方向上的位移；Δz(x,y)為變動平面z(x,y)上任意兩點的z坐標值之差。由于變動后公差實際表面的形狀保持不變,仍為矩形平面,極值情況必然發生在平面的4個頂點處,因此只需考慮變動平面頂點的z向坐標值是否在公差域內即可。

由圖1可知,4個頂點的x、y向坐標分別為(a,b)、(-a,b)、(a,-b)、(-a,-b),由實際變動平面方程可知,4個頂點的坐標分別為

S1=(a,b,dz+aβ+bα)

S2=(-a,b,dz-aβ+bα)

S3=(a,-b,dz+aβ-bα)

S4=(-a,-b,dz-aβ-bα)

由式(1)可得

(3)

由式(2)、式(3)可得旋量參數的變動不等式[8]為

(4)

約束不等式[8]為

(5)

-TL≤w+xβ+yα≤TU

(6)

其中，x、y分別取4個頂點的x、y向坐標值。

變動不等式式(4)定義了旋量參數的理想變動區間。約束不等式的存在使旋量參數不可能同時取最大值，參數的實際變動區間比理想變動區間要小，如果按旋量參數的理想變動區間進行精度設計,則會提高零件加工精度的要求，增加制造成本。而旋量參數實際變動區間與公差間的具體關系尚不明確。

2.1.2蒙特卡洛模擬法求解旋量參數實際變動區間

蒙特卡洛模擬法也稱為隨機模擬方法,它以概率論和數理統計為基礎，通過對隨機變量的統計實驗和隨機模擬來求解問題近似解。由于多種因素綜合作用的影響，產品公差實際表面變動的形式在公差域內是無法預知的，具有隨機性，其分布的規律隨加工條件而定,因此,與公差對應的SDT旋量參數也是隨機量[13]。由式(4)～式(6)可知,旋量參數的變動順序不同，其最終的取值也會不同。本文根據變動不等式和約束不等式，運用蒙特卡洛模擬法，按給定的參數分布類型和變動順序進行模擬試驗，生成相應的參數隨機數，模擬實際公差表面的變動，對滿足約束條件的隨機數加以保留，不滿足約束條件的則剔除，當產生足夠多滿足條件的隨機數樣本后，對樣本進行分析，求解出對應旋量參數的變動區間帶寬，再根據變動區間帶寬和均值得到旋量參數的實際變動區間。求解步驟如下:

(1)根據研究對象的具體情況，確定旋量參數的理想概率分布模型。本文假定旋量參數的分布函數均為正態分布:

式中，μ、σ分別為正態分布的均值和均方差。

(3)根據抽樣規則對旋量參數進行抽樣。3個非零旋量參數共有6種變動順序，為了盡可能地模擬零件的實際加工工況,對每一種變動順序都進行抽樣。以變動順序α→β→w為例,其抽樣流程如圖2所示。

圖2　抽樣流程一

圖2中,k的初始值為1;K為要求抽取的合格樣本數,為保證模擬精度,K>20 000;K1、K2為常量,根據具體研究對象而定。設置判別式k1

按6種變動順序編程進行抽樣后,α、β、w均得到容量為6K的實際變動區間帶寬樣本。

(4)旋量參數實際分布函數的假設檢驗。由于約束不等式式(5)、式(6)的限制，旋量參數的實際分布類型并不一定與理想分布類型相同,故需對實際分布類型進行假設檢驗。本文采用χ3擬合檢驗法,以參數α為例,其基本思路為:在α分布未知時,先根據抽樣獲得的觀測值對α的分布類型作出假設:

H0∶F(α)=F0(α)

式中,F0(α)為假設的α分布函數。

將實數軸分為k個互不相交的區間(ai,ai+1](i=1,2,…,k),其中a1、ak+1可分別取-∞、+∞,區間的劃分視具體情況而定。在假設H0下計算概率:

pi=F(ai+1)-F(ai)=P(ai<α

計算理論頻數。由于α的樣本容量為6K,因此理論頻數為6Kpi。

計算樣本觀察值落入區間(ai,ai+1]中的個數fi(i=1,2,…,k),fi稱為實際頻數。

計算統計量

當樣本容量6K≥50時,統計量χ2近似服從自由度為k-r-1的χ2分布(r為確定F0(α)表達式中被估計參數的個數,如F0(α)為正態分布,則r=2)。

(5)采用極大似然估計法估計實際分布函數的參數。如果經第(4)步假設檢驗驗證α也屬于正態分布N(μ,σ2),則根據α的樣本,采用極大似然估計法估計α分布函數的參數,得到μ與σ2的極大似然估計量為

由于變動不等式與約束不等式均具有對稱性,故旋量參數實際分布的均值近似等于理想分布的均值。

(6)求解旋量參數的實際變動區間。根據參數的實際分布函數類型與均方差，查表得到參數的實際變動區間帶寬:

其中,G表示相對分布系數,如果旋量參數為正態分布，則G等于1。旋量參數的實際變動區間為

2.1.3運用響應面法建立旋量參數實際變動區間帶寬與公差間的響應面函數

響應面方法是數學方法與統計方法相結合的一種建模方法,主要用于解決如何建立復雜系統輸入(變量)與輸出(響應)之間近似函數關系的問題。其基本思路是,依據試驗設計原則選定設計空間中的試驗點，用多組試驗點及其對應的響應值擬合響應面曲面，從而構造具有明確表達形式的顯式函數[14-15]。

通過前文分析可知，只要確定了旋量參數實際變動區間帶寬即可確定參數的實際變動區間。對于任意的平面尺寸公差均可通過蒙特卡洛模擬法得到與之對應的旋量參數實際變動區間帶寬的響應值，但兩者間的函數關系并不明確，無法用嚴格的數學公式表達。為此，本文采用響應面方法，建立公差與旋量參數實際變動區間帶寬間的響應面模型,進而確定公差與旋量參數實際變動區間之間的響應關系，為后續的精度分析與精度分配研究打下基礎。建模的步驟如下:

Fy=c0+c1T+c2T2y=α,β,w

式中,ci(i=0,1,2)為待求系數。

采用最小二乘法求解C=[c0c1c2]T的無偏估計:

(3)驗證模型精度。響應面生成后，采用復相關系數R2對響應面模型進行驗證，其表達式為

復相關系數R2代表了響應面預測值與真值之間的差異程度,在0～1之間取值,值越大,表示差異度越小。

2.2平面度公差建模

形位公差數學建模的目的是確定形位公差域邊界的變動范圍及公差變動要素的表示方法。由工程實踐可知，對某幾何要素給出形位公差的原因是：當其他的公差(如尺寸公差)存在，但仍不能滿足對幾何要素形位變動限制的要求時，需要給出公差值較小的形位公差來進一步限制。實際上，尺寸公差與形位公差是同一實際幾何要素反映的兩種不同概念的公差，是同一問題的兩個方面。因此，在研究形位公差的建模時不能離開尺寸公差而單獨處理，需要進一步研究形位公差與尺寸公差的關系，即要考慮公差原則。

由GB4249-1996可知，公差原則包含獨立原則和相關要求，相關要求中又分為包容要求、最大實體要求和最小實體要求。本文以獨立原則為例，對平面度公差進行建模。

2.2.1求解平面度公差變動不等式與約束不等式

由平面度的定義可知,若平面度大小為TF,則平面上所有的點均必須位于距離為TF的兩平行平面所形成的公差域內。用Bottom和Top分別表示兩平行平面,Bottom面表示材料邊內的平面，Top面表示材料外的平面,平面度公差規定的形式如圖3所示。

(a)平面度公差

(b)平面度公差域圖3　平面度公差及其公差域

圖3中,+TSU、-TSL分別表示尺寸公差TS的上下偏差;TF表示平面度公差;平面A為基準平面,平面B為變動平面。設Bottom與Top平面對應的SDT分別為(αB,βB,0,0,0,wB)T、(αT,βT,0,0,0,wT)T,則Bottom面和Top面的變動方程分別為

zB=wB+xβB+yαB

zT=wT+xβT+yαT=

由獨立原則可知,平面B必須在尺寸公差域(ST域)內,即與平面A相距D-TSL～D+TSU的區域。此時,對平面度公差域(FT域)的唯一限制是在FT域內存在一個平面，它同時也在ST域。采用基于數學定義的公差分析方法求解得到Bottom面旋量參數的變動不等式為

(7)

約束不等式為

-TSL-TF≤xβB+yαB≤TSU

(8)

-TSL-TF≤wB+xβB+yαB≤TSU

(9)

其中x、y分別取Bottom面4個頂點的x、y向坐標值。

由于Bottom面和Top面互相平行,故有

αT=αB,βT=βB,wT=wB+TF

(10)

通過式(7)～式(10)確定了Bottom和Top面的變動范圍,但變動平面在Bottom和Top面限制的區域中如何分布并不確定。為此,本文采用二元線性回歸法求解變動平面旋量參數的值,其基本思路如下:已知Bottom和Top面的方程，在名義平面內設置m×n個取樣點，使之成m×n網格分布，通過隨機模擬法對每個取樣點在Bottom和Top面限制的范圍內取z值,這樣就得到Q(Q=m×n)個平面度區域內的模擬點。根據Q個模擬點的值，采用二元線性回歸算法計算得到變動平面旋量參數的值。旋量參數的表達式如下:

2.2.2蒙特卡洛模擬法求解旋量參數實際變動區間

由于可以假設變動后公差實際表面的形狀保持不變,故用蒙特卡洛模擬法求解尺寸公差旋量參數實際變動區間時，通過模擬抽樣得到的旋量參數值即可作為公差的實際旋量參數。與尺寸公差不同之處在于，形狀公差只確定了公差特征要素的分布區域，通過蒙特卡洛模擬法抽樣只能確定公差特征要素分布區域的上下邊界(即Bottom和Top面),而特征要素在區域內如何分布仍不確定，需要再采用二元線性回歸法求解公差的實際旋量參數。求解步驟如下：

(1)根據研究對象的具體情況，確定旋量參數的理想概率分布模型。本文假定旋量參數的分布函數均為正態分布:

(2)根據Bottom面變動不等式式(7)確定各旋量參數理想分布的均值與方差。由式(7)可知，旋量參數αB、βB、wB理想分布的均值與方差分別為

(3)根據抽樣規則對旋量參數進行抽樣。Bottom面的3個非零旋量參數同樣有6種變動順序,需對每一種變動順序都進行抽樣。以變動順序αB→βB→wB為例,其抽樣流程如圖4所示。

圖4　抽樣流程二

按6種變動順序進行抽樣后,α、β、w均得到容量為6K的實際變動區間帶寬樣本。

第(4)～(6)步與2.1.2節所述過程相似，在此不再重復闡述。

2.2.3運用響應面法建立旋量參數實際變動區間帶寬與公差間的響應面函數

由式(7)可知,與尺寸公差響應面方法建模的不同之處在于，平面度公差建模過程中包含尺寸公差和平面度公差2個設計變量,需同時考慮這兩者對變動區間帶寬的影響,建模的步驟如下：

式中,ci(i=0,1,…,4)為待求系數。

采用最小二乘法求解C=[c0c1c2c3c4]T的無偏估計

(3)驗證模型精度。與2.1.3節所述過程相似,在此不再重復闡述。

3　實例分析

圖5　α變動區間帶寬

圖6　β變動區間帶寬

圖7　w變動區間帶寬

令下降比例等于殘余帶寬除以理想變動區間帶寬。下降比例均值為100個試驗點對應的下降比例的均值。下降比例均值反映了將蒙特卡洛模擬法應用于精度設計時對提高精度設計技術經濟性的效果。表2列出了3個旋量參數實際變動區間帶寬下降比例均值。

表2　實際變動區間帶寬下降比例均值　%

由表2可知,在考慮約束條件及參數變動順序后，運用蒙特卡洛模擬法求解得到的旋量參數實際變動區間帶寬較理想變動區間帶寬下降比例的均值分別為22.49%、22.85%、20.34%,說明采用蒙特卡洛分析方法可以顯著提高精度設計的技術經濟性,且更加符合零件生產實際。

Fα=0.0026+0.024TS+

Fβ=0.0014+0.188TS+

Fw=0.0657+0.57TS+

用復相關系數R2驗證3個響應面模型的精度,如表3所示。由表3可知,3個響應面模型的R2值均大于0.985，表明本文建立的公差響應面

表3　響應面模型的R2檢驗

模型具有很高的精度。

4　結語

本文將蒙特卡洛模擬法良好的求解非確定性問題的能力與響應面方法優秀的建模性能結合起來，提出了一種基于蒙特卡洛模擬與響應面方法的公差建模方法。在對平面尺寸和形狀公差域變動關系進行深入研究的基礎上，采用蒙特卡洛模擬法求解公差與公差變動要素變動區間之間的響應關系，再運用響應面方法建立了兩者之間的響應面數學模型。以平面度公差建模為例，驗證了該方法可以顯著提高精度設計的技術經濟性。進一步將該方法應用到機床的精度設計當中，對提高機床精度設計方法的實用性,降低機床的生產成本具有重要意義。

[1]劉又午.多體動力學的休斯敦方法及其發展[J].中國機械工程,2000,11(6):601-608.

Liu Youwu.Development of Huston’s Method on Multibody Dynamics[J].China Mechanical Engineering,2000,11(6):601-608.

[2]Chen J S.Computer-aided Accuracy Enhancement for Multi-axis CNC Machine Tool[J].International Journal of Machine Tools and Manufacturing,1999,35(4):593-605.

[3]Patel A J,Ehmann K F.Volumetric Error Analysis of a Stewart Platform-based Machine Tool[J].Annals of the CIRP,1997,46(1):287-290.

[4]Laperriere L,Ghie W,Desrochers A.Statistical and Deterministic Tolerance Analysis and Synthesis Using a Unified Jacobian-Torsor Model[J].CIRP Annals-Manufacturing Technology,2002,51(1):417-420.

[5]張為民，李國偉，陳燦．基于雅可比旋量理論的公差優化分配[J]．農業機械學報，2011，42(4)：216-219.

Zhang Weimin,Li Guowei,Chen Can.Optimal Allcocation of Tolerance Based on Jacobian-Torsor Theory[J].Transactions of the Chinese Society for Agricultural Machinery,2011,42(4)：216-219.

[6]胡潔,吳昭同,楊將新.基于旋量參數的三維公差累積的運動學模型[J].中國機械工程,2003,14(2):127-130.

Hu Jie,Wu Zhaotong,Yang Jiangxin.Kinematic Model of 3D Tolerance Accumulation Based on Screw Parameter[J].China Mechanical Engineering,2003,14(2):127-130.

[7]Mujezinovic A, Davidson J K, Shah J J. A New Mathematical Model for Geometric Tolerances as Applied to Polygonal Faces[J]. Transections of the ASME,2004,126(3):504-518．

[8]劉玉生,吳昭同,楊將新,等.基于數學定義的平面尺寸公差數學模型[J].機械工程學報,2001,37(9):12-17.

Liu Yusheng,Wu Zhaotong,Yang Jiangxin,et al.Mathematical Model of Size Tolerance for Plane Based on Mathematical Definition[J].Chinese Journal of Mechanical Engineering,2001,37(9):12-17.

[9]Roy U,Li B.Representation and Interpretation of Geometric Tolerances for Polyhedral Objects Form Tolerances[J].Computer Aided Design,1998,30(2):151-161.

[10]蔡敏,楊將新,吳昭同.平面要素公差數學定義理論及應用的研究[J].中國機械工程,2002,13(2):128-130.

Cai Min,Yang Jiangxin,Wu Zhaotong.Theory and Application of Mathematical Definition for Plane Element Tolerance[J].China Mechanical Engineering,2002,13(2):128-130.

[11]Vignat F,Villeneuv F.3D Transfer of Tolerances Using a SDT Approach:Application to Turning Process[J].Journal of Computing and Information Science in Engineering,2003,3(3):45- 53.

[12]Shah J J.Dimension and Tolerance Modeling and Transformations in Feature Based Design and Manufacturing[J].Journal of Intelligent Manufacturing,1998(9):475-488．

[13]Zhong Xin,Yang Ruqing,Zhou Bing.Accuracy Analysis of Assembly Success Rate with Monte Carlo Simulations[J].Journal of Donghua University,2003,20(3):128-131.

[14]Kim S H,Na S W.Response Surface Method Using Vector Projected Sampling Points[J].Structrual Safety,1997,19:3-19.

[15]Kaymaz I,Mcmahon A.A Response Surface Method Based on Weighted Regression for Structural Reliability Anslysis[J].Probabilistic Engineering Mechanics,2004,20：11-17.

(編輯蘇衛國)

Tolerance Modeling Based on Monte-Carlo Simulation and Response Surface Method

Yu ZhiminLiu ZijianDong SikeLi SiminAi Yandi

State Key Laboratory of Advanced Design and Manufacture for Vehicle Body,Hunan University,Changsha,410082

For the poor practicability of existing precision design method of machine tool,a new method about tolerance modeling based on MCS and response surface method was proposed herein.Firstly,the constraint inequality and change inequality were established using tolerance analysis based on mathematic definition theory.And then the simulated tests with MCS method were conducted to simulate the change of actual tolerance surface,and get the actual change interval of the tolerance factor.With the samples based on the tolerance and the bandwidth value of actual change interval obtained above,the response surface modeling between the tolerance and bandwidth was constructed by response surface method.Finally,a typical example was analyzed and the results indicate the method is consistent with engineering practice, and has high modeling precision and technical economics.

Monte-Carlo simulation(MCS);response surface method;small displacement torsor;constraint inequality;change inequality

2013-06-21

國家自然科學基金資助項目(51175161,51301532);國家科技重大專項(2011ZX04003-011)

TH115DOI：10.3969/j.issn.1004-132X.2015.04.001

余治民,男,1984年生。湖南大學汽車車身先進設計制造國家重點實驗室博士研究生。研究方向為復雜裝備精度鏈設計方法。劉子建,男,1953年生。湖南大學汽車車身先進設計制造國家重點實驗室教授、博士研究生導師。董思科,男,1986年生。湖南大學汽車車身先進設計制造國家重點實驗室碩士研究生。李斯明，男，1988年生。湖南大學汽車車身先進設計制造國家重點實驗室碩士研究生。艾彥迪,男,1982年生。湖南大學汽車車身先進設計制造國家重點實驗室助理研究員。