基于PCE的諧波減速器動態精度不確定性分析

張金洋， 張建國，*， 彭文勝， 劉育強， 汪龍

(1. 北京航空航天大學 可靠性與系統工程學院， 北京 100083； 2. 北京航空航天大學 可靠性與環境工程重點實驗室， 北京 100083； 3. 中國航空綜合技術研究所， 北京 100028； 4. 中國空間技術研究院 總體部， 北京 100029；5. 北京衛星制造廠， 北京 100094)

諧波減速器是一種建立在彈性變形理論上的新型機械傳動方式，在航空航天、機器人等精密定位領域廣泛應用。相比于其他傳統的傳動方式，諧波傳動具有精度高、傳動比大、傳動平穩和傳動效率高等優點。為了提高諧波減速器的傳動精度，國內外學者對諧波傳動誤差進行了深入的研究[1-4]，主要從諧波減速器的加工、裝配因素等傳動機理方面考慮，并沒有考慮諧波減速器的柔性和非線性摩擦等動態特性對其傳動誤差的影響。另外，在對諧波減速器進行靜態誤差分析時，并沒有考慮相關參數不確定性的影響以及靜態誤差和動態誤差之間的耦合關系。Hsia[5]提出了諧波傳動誤差主要是由波發生器帶動柔輪運動時柔輪變形引起的，從設計的角度研究柔輪柔性變形對動態誤差的影響，但未考慮其他非線性因素的影響。Tuttle和Seering[6]考慮了諧波減速器柔性以及靜態誤差等非線性因素，同時深入研究了剛柔輪齒嚙合的非線性機理，建立了非線性動力學模型，并通過試驗驗證了所建模型的正確性。游斌弟和趙陽[7]研究了諧波減速器動態誤差在考慮齒輪嚙合摩擦和扭轉剛度非線性因素時的動力學響應，利用拉格朗日方程建立諧波減速器的誤差動力學方程，研究了不同頻率諧波強迫激勵下動態誤差影響。Preissner等[8]建立了綜合考慮柔輪非線性扭轉、滯回特性和運動誤差的諧波傳動模型，著重研究了柔輪的滯回特性對動態精度的影響。這些研究大多只考慮諧波減速器的柔輪動力學特性對動態精度的影響，沒有考慮動力學模型參數的不確定性。

本文從靜態誤差和動態誤差產生的機理方面著手，綜合考慮了加工、裝配誤差、柔輪柔性以及剛柔輪摩擦綜合作用下的諧波減速器的動態精度問題；在此基礎上考慮相關參數的不確定性對動態精度的影響，利用多項式混沌展開(Polynomial Chaos Expansion，PCE)方法進行動態精度不確定性分析，構建動態精度可靠度模型計算動態精度可靠度。結果表明，利用PCE方法能夠很好地處理諧波減速器非線性以及多變量的動態精度不確定性問題。

1 諧波減速器動態精度分析

諧波減速器主要由3部分組成：波發生器、柔輪和剛輪，其結構簡圖如圖1 所示。電機輸出一定轉速驅動波發生器旋轉，波發生器帶動柔輪發生柔性變形與剛輪嚙合傳遞運動。其中波發生器是輸入端，柔輪是輸出端。

(1)

式中：θm為波發生器輸入角位置；θl為柔輪輸出角位置；N為諧波傳動比。

在一定轉速下諧波減速器柔輪會發生柔性變形，當其做正反往復運動時，就會出現如圖2 所示的滯回誤差。

圖1 諧波減速器結構簡圖Fig.1 Structure diagram of harmonic reducer

圖2 諧波減速器往復運動動態誤差曲線Fig.2 Dynamic error curve of harmonic reducer reciprocation

(2)

靜態誤差[10]主要由諧波減速器各個部件的加工、裝配誤差所產生。其中剛輪加工誤差產生的運動誤差Δc1為

(3)

柔輪加工誤差引起的運動誤差Δr1為

(4)

剛輪、柔輪裝配誤差，不考慮裝配誤差初相角的影響，由剛輪安裝偏心誤差Ec產生的運動誤差為

Δc2=Ecsin(2ωbt)/cosαn

(5)

式中：αn為嚙合角。由柔輪安裝偏心誤差Ef產生的運動誤差為

(6)

由波發生器安裝偏心誤差Eb產生的運動誤差Δb為

Δb=Ebsin(ωbt)/cosαn

(7)

綜上，可以得到總的靜態誤差為

(8)

式中：d為剛輪分度圓直徑。

動態精度包括機構的靜態誤差和柔輪柔性和機構摩擦引起的滯回誤差，分析諧波減速器的動態精度，應該在靜態誤差模型基礎之上建立機構的動力學模型。本文研究的諧波減速器物理簡化模型如圖3所示，圖中：Jm和Jl分別為輸入和輸出端轉動慣量；Bm、Bl和Bsp分別為諧波減速器輸入端、輸出端阻尼和剛柔嚙合處阻尼。

(9)

圖3 諧波減速器物理簡化模型Fig.3 A physical simplified model of harmonic reducer

式中：Tk為柔輪的非線性扭矩；k1和k2為諧波減速器柔輪等效扭轉剛度。由圖3所示的諧波減速器各部件之間的動力學關系，系統的動能T為

(10)

系統的勢能V為

(11)

系統的瑞利耗散函數為

(12)

因此，諧波減速器拉格朗日動力學方程[11]為

(13)

(14)

(15)

式中：Tm為輸入轉矩。

2 基于PCE動態精度不確定性模型

2.1 PCE基本理論

PCE基本理論是用一個屬于某個對應分析類型的正交多項式混沌之和(含有一個或多個隨機變量)來近似地表示一個隨機過程。

(16)

c=(c0,c1,…)為待定系數矢量；ξ=[ξ1,ξ2,…,ξn]為服從標準正態分布的隨機變量矢量；Πd(ξi1,ξi2,…,ξid) 為d次多維Hermite多項式。將式(16)截斷用s項來近似精度則可簡化[12]為

(17)

式中：η為隨機事件；cj為待求解的確定性系數；Πj(ξ1,ξ2,…,ξn)為廣義Wiener-Askey多項式混沌。對于一個n維Hermite多項式，可表示為

(18)

Hermite多項式的隨機變量如果是標準正態分布，則滿足

〈Πp(ξ)Πq(ξ)〉=0p≠q

(19)

式中：〈··〉 為希爾伯特空間上的內積，此處定義為

(20)

式中：多項式基對應的權重函數W(ξ)為

(21)

如果機構輸入的隨機變量的個數為n，又機構輸出響應量的多項式展開式最高階次為p，則待定系數的個數P可以用式(22)求得：

(22)

2.2 動態精度不確定性模型

諧波減速器的動態精度不僅跟機構部件的制造公差、裝配間隙等參數有關，而且還受機構動力學參數(如等效剛度、阻尼和轉動慣量)的影響。上述參數在制造、測量過程中必然會存在不確定性而不是一個固定的數值，參數的不確定性通過機構的動力學模型影響動態精度響應的不確定性。

設諧波減速器的動力學模型用M表示，機構部件的制造公差、間隙、等效剛度和轉動慣量等不確定參數可以用向量ψ表示，即

ψ=[ψ1,ψ2,…,ψn]

(23)

(24)

(25)

(26)

(27)

式中：ci,m、cij,m、cijj,m和cijk,m為展開式中的待定系數。

本文采用隨機響應面配點法[14]計算多項式混沌的展開系數。在隨機向量展成的空間中，每一組樣本{ξ1,ξ2，…}都會對應一個點，這些點稱為配點，對于Wiener的多項式混沌，其展開式基函數為Hermite多項式 ，若Hermite多項式的最高階數為p，則相應的配點通常取p+1階Hermite多項式的根。未知系數可以通過式(28)和式(29)計算：

[c0(t)c1(t) …cP-1(t)]T=(Π(ξ)TΠ(ξ))-1·

Π(ξ)T[θ(t,ξ0)θ(t,ξ1) …θ(t,ξk)]T

(28)

(29)

式中：ξ1,ξ2,…,ξk為采樣點；k為采樣點數。

根據Hermite多項式正交性，動態精度的均值可以通過式(30)求得[15]：

(30)

(31)

3 不確定性及可靠性求解分析

3.1 諧波減速器Dymola仿真模型

根據諧波減速器動力學模型，本文針對XB1-50型號諧波減速器采用Modelica語言在Dymola[16]編譯環境下建立諧波減速器的動力學仿真模型(見圖4),主要包括直PID控制模塊、電機模塊和諧波減速器模塊。此模型能夠較好處理諧波減速器的機、電耦合問題，在求解非線性微分方程時精度較高。

采用如圖5所示的諧波減速器精度測試平臺對所建立的非線性諧波減速動力學仿真模型進行驗證。采用恒溫箱封閉，電機輸入端和諧波減速器輸出端分別連接有角速度和角度傳感器，測定輸入、輸出端角位置以及角速度。實驗條件下，電機輸入轉速為100 r/min。

根據實測和設計參數數據確定諧波減速器動力學模型各參數如表1所示，靜態誤差模型參數如表2所示。

在Dymola仿真模型中通過PID調節器控制電機輸出轉速為100 r/min，仿真時間為10 s，待電機輸入轉速及誤差波動曲線穩定后，通過采樣輸入角在3個周期均勻變化下動態誤差值，實驗條件下，通過角位移傳感測得同樣周期內諧波減速器輸入端角位移θm以及輸出端角位移θl，求得動態誤差的測量值。將動態誤差仿真值與測量值對比如表3所示。

圖4 諧波減速器Dymola仿真模型Fig.4 Dymola simulation model of harmonic reducer

圖5 諧波減速器精度測試平臺Fig.5 Precision test platform of harmonic reducer

參 數數 值Jm/(kg·m2)3.2×10-4Jl/(kg·m2)8.5×10-4Bm/(N·m·s·rad-1)1.7×10-4Bl/(N·m·s·rad-1)5.0×10-4Bsp/(N·m·s·rad-1)2.8×10-4k1/(N·m·rad-1)7160k2/(N·m·rad-3)21576N90

表2 靜態誤差模型參數

表3 不同電機輸入角下動態誤差實驗值與仿真值

通過表3中數據可以看出，實驗得到的動態誤差值與仿真值比較接近，相對誤差在2.5%～9.5%范圍內，考慮到相關參數不確定性，相對誤差在接受范圍之內，由此可見所建立的諧波減速器Dymola仿真模型能夠很好地模擬考慮間隙和柔性非線性因素綜合作用下的動態精度。

仿真得到諧波減速器輸入、輸出轉速，動態誤差和靜態誤差曲線如圖6～圖9所示， 圖6為諧波減速器輸入轉速曲線，轉速經過短暫波動后穩定后在10.499 rad/s。由于諧波減速器柔性和控制器慣性的作用，圖7所示諧波減速器輸出轉速先出現大范圍波動，然后穩定在0.1 rad/s附近波動。圖8為動態誤差隨輸入角θm波動曲線，開始波動程度較大，穩定后，呈周期性波動。由圖9可以看出，綜合考慮靜態誤差、柔性和摩擦作用的動態誤差曲線比只考慮靜態誤差曲線波動幅值要大0.005°，相比增加25.4%，因此在對諧波減速器進行動態誤差分析時，有必要考慮柔性和摩擦的影響。

圖6 諧波減速器輸入轉速Fig.6 Speed of harmonic reducer input

圖7 諧波減速器輸出轉速Fig.7 Speed of harmonic reducer output

圖8 諧波減速器動態誤差Fig.8 Dynamic error of harmonic reducer

圖9 靜態誤差和動態誤差曲線Fig.9 Static error and dynamic error curves

3.2 動態精度靈敏度及不確定性分析

Sobol敏感度分析方法[17]是一種基于方差分解的Monte Carlo方法，Sobol方法考慮了隨機輸入參量在整個取值范圍內對輸出響應的貢獻以及隨機參數的交互作用，Sobol敏感度指標可以直接從PCE式的系數得到。

動態精度PCE式(26)和式(27)可以改寫為

(32)

式中:cα為PCE多項式系數；ψα為PCE多項式的項；Ii1i2…is={α∈(α1,α2,…,αN):αk=0?k?(i1i2…is),?k=1,2,…,N},根據多項式混沌的正交性，可得

(33)

(34)

由式(34)可以看出分解項表征不同隨機變量及其相互作用對諧波減速器動態精度輸出響應方差的貢獻，因此可以定義Sobol敏感度指標為

(35)

其滿足

(36)

式中：Si為主效應敏感度指標，表征各個隨機變量對動態精度響應方差的貢獻。

表4 諧波減速器不確定性參數及分布

表5 動態精度多項式混沌展開式配點

表6 動態精度多項式混沌展開式系數

圖10 不同扭轉剛度k1時動態誤差變化曲線Fig.10 Dynamic error curves at different torsional stiffness k1

圖11 不同輸出軸轉動慣量Jl時動態誤差變化曲線Fig.11 Dynamic error curves at different output shaft moment of inertia Jl

圖12 不同柔輪切向相鄰齒綜合誤差時動態誤差變化曲線Fig.12 Dynamic error curves at different flexible wheel tangential adjacent gear comprehensive error

圖13 2種方法動態誤差均值比較Fig.13 Dynamic error mean comparison between two methods

3.3 動態精度可靠性分析

本文定義諧波減速器動態精度可靠度為諧波減速器在一定初始條件和時域內，動態誤差絕對值的最大值不超過一定的閾值的概率，則動態精度可靠性功能函數可表示為

(37)

圖14 2種方法動態誤差均方差比較Fig.14 Dynamic error’s mean square error comparison between two methods

項系 數項系 數項系 數10.024151ζ2-1-0.0008ζ1ζ5-0.0053ζ1-0.2361ζ3-10.0490ζ2ζ3-0.0054ζ20.0366ζ4-1-0.0002ζ2ζ40.0001ζ30.7209ζ5-1-0.0015ζ2ζ50.0015ζ4-0.0236ζ1ζ20.0037ζ3ζ40.0039ζ5-0.0218ζ1ζ3-0.0365ζ3ζ50.0042ζ1-10.0365ζ1ζ40.0006ζ4ζ5-0.0005

(38)

(39)

通過迭代計算得到動態精度可靠度計算結果如表8所示。

在Dymola中進行500次Monte Carlo抽樣仿真實驗，得到動態誤差仿真曲線如圖15所示。

根據式Pf=nf/nt求動態精度在一定時域內的失效概率，nf為動態誤差絕對值最大值超過閾值的次數，nt為抽樣仿真次數，表9給出了nt=500～10 000次的模擬結果。

從表9中可以看出,當n=8 000時，Pf收斂，Monte Carlo仿真實驗得到動態精度可靠度為0.961 9，利用PCE方法和FORM法得到的可靠度0.958 4，兩者求得的結果相近，而PCE方法僅用到42次抽樣仿真，效率明顯優于Monte Carlo方法。

表8 FORM法可靠度計算結果

圖15 動態誤差的Monte Carlo仿真曲線Fig.15 Monte Carlo simulation curves of dynamic error

ntnfPf500140.02802000610.030540001260.031560002050.034270002540.036380003050.0381100003810.0381

4 結 論

本文考慮了靜態誤差、柔性和摩擦等因素綜合作用下的諧波減速器動態精度問題，建立了諧波減速器非線性動力學模型，利用多項式混沌展開(PCE)方法處理動態精度的隨機不確定性。解決了以下問題：

1) 將靜態誤差與動態誤差結合起來，建立了含有靜態誤差項的非線性動力學模型，通過實驗驗證了考慮動力學因素和靜態誤差綜合作用下的動態精度更接近實際情況。

2) 基于PCE方法，通過靈敏度分析得到影響動態精度的主要參數是柔輪等效剛度、輸出軸轉動慣量和柔輪切向相鄰齒綜合誤差?？梢钥闯?，柔輪的加工精度對動態精度影響較大。通過靈敏度分析選取影響較大的參數進行不確定性分析，可以減少不確定性分析的計算量。

3) 利用PCE方法得到動態精度的隨機統計特性，并將PCE方法與傳統的Monte Carlo方法比較，效率更高。利用PCE方法和FORM法計算一定時域內的動態精度可靠度，比傳統的Monte Carlo方法計算可靠度效率更高。

參考文獻 (References)

[1] NYE T,KRAML R.Harmonic drive gear error:Characterization and compensation for precision pointing and tracking[C]∥Processing of the 25th Aerospace Mechanics,Symposium.Washington,D.C.：NASA，1991:237-252.

[2] EMELYANOV A F.Calculation of the kinematic error of a harmonic gear transmission taking into account the compliance of elements[J].Soviet Engineering Research,1983,3(7):7-10.

[3] GRAVAGNO F,MUCINO V H,PENNESSTRI E.Influence of wave generator profile on pure kinematic error and centrodes of harmonic drive[J].Mechanism and Machine Theory,2016,104:100-107.

[4] 沙曉晨，范元勛.諧波減速器傳動誤差的研究[J].機械制造,2015,44(5):50-54.

SHA X C，FAN Y X.Study of transmission error of harmonic drive reducer[J].Machine Building Automation,2015,44(5):50-54(in Chinese).

[5] HSIA L M.The analysis and design of harmonic gear drives[C]∥Processing of the 1988 IEEE International Conference on Systems,Man,and Cybernetics.Piscataway,NJ:IEEE Press,1988:616-620.

[6] TUTTLE T,SEERING W.Kinematic error,compliance,and friction in a harmonic drive gear transmission[C]∥ASME Design Technical Conferences 19th Design Automation.New York:ASME,1993:319-324.

[7] 游斌弟，趙陽.考慮非線性因素的諧波齒輪傳動動態誤差研究[J].宇航學報，2010，31(5):1297-1282.

YOU B D,ZHAO Y.Study on dynamic error of harmonic drive with nonlinear factors[J].Journal of Astronautics,2010,31(5):1297-1282(in Chinese).

[8] PREISSNER C,ROYSTON T J,SHU D.A high-fidelity harmonic drive model[J].Journal of Dynamic Systems,Measurement,and Control,2012,134(1):011002.

[9] KIRCANSKI N,GOLDENBERG A,ANGELS J.Nonlinear modeling and parameter identification of harmonic drive gear transmissions[C]∥Processing of the 32nd IEEE Conference on Robotics and Automatic.Piscataway,NJ:IEEE Press,1995:3027-3032.

[10] 王愛東.機器人用諧波齒輪傳動裝置的運動精度分析[D].北京:中國科學院,2001:11-15.

WANG A D.Kinematic accuracy analysis of gear transmission for harmonic reducer of robot[D].Beijing:Chinese Academy of Sciences,2001:11-15(in Chinese).

[11] FATHI H,PRASANNA S,FRIEDHELM A.On the kinematic error in harmonic drive gears[J].Journal of Mechanical Design,2001,123(1):90-97.

[12] WIENER N.The homogeneous chaos[J].American Journal of Mathematics,1938,60(4):897-936.

[13] 趙珂,高正紅,黃江濤,等.基于PCE方法的翼型不確定性分析及穩健設計[J].力學學報,2014,46(1):11-19.

ZHAO K,GAO Z H,HUANG J T,et al.Uncertainty quantification and robust design of airfoil based on polynomial chaos technique[J].Chinese Journal of Theoretical and Applied Mechanics,2014,46(1):11-19(in Chinese).

[14] LOEVETT T,PONCI F,MONTI A.A polynomial chaos approach to measurement uncertainty[J].IEEE Transations on Instrumentation and Measurement,2006,55(3):729-736.

[15] BENJAMIN L,HOSAM K,JEFFREY L.Recursive maximum likelihood parameter estimation for state space systems using polynomial chaos theory[J].Automatica,2011,47(11):2420-2424.

[16] 陶海川,來新民.基于Dymola 的無刷直流電機仿真模型[J].計算機仿真，2005,22(5):63-65.

TAO H C,LAI X M. Computer simulation of brushless DC motor system based on Dymola[J].Computer Simulation，2005,22 (5):63-65(in Chinese).

[17] SUDRET B.Global sensitivity analysis using polynomial chaos expansion[J].Reliability Engineering and System Safety,2008,93(7):964-979.

[18] PENG W S,ZHANG J G,ZHU D T.ABCLS methods for high-reliability aerospace mechanism with truncated random uncertainties[J].Chinese Journal of Aeronautics,2015,28(4):1066-1075.