基于PCE的諧波減速器動態(tài)精度不確定性分析

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

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

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

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

1 諧波減速器動態(tài)精度分析

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

(1)

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

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

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

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

(2)

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

(3)

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

(4)

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

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

(5)

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

(6)

由波發(fā)生器安裝偏心誤差Eb產(chǎn)生的運動誤差Δb為

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

(7)

綜上，可以得到總的靜態(tài)誤差為

(8)

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

動態(tài)精度包括機構的靜態(tài)誤差和柔輪柔性和機構摩擦引起的滯回誤差，分析諧波減速器的動態(tài)精度，應該在靜態(tài)誤差模型基礎之上建立機構的動力學模型。本文研究的諧波減速器物理簡化模型如圖3所示，圖中：Jm和Jl分別為輸入和輸出端轉(zhuǎn)動慣量；Bm、Bl和Bsp分別為諧波減速器輸入端、輸出端阻尼和剛?cè)釃Ш咸幾枘帷?/p>

(9)

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

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

(10)

系統(tǒng)的勢能V為

(11)

系統(tǒng)的瑞利耗散函數(shù)為

(12)

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

(13)

(14)

(15)

式中：Tm為輸入轉(zhuǎn)矩。

2 基于PCE動態(tài)精度不確定性模型

2.1 PCE基本理論

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

(16)

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

(17)

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

(18)

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

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

(19)

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

(20)

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

(21)

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

(22)

2.2 動態(tài)精度不確定性模型

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

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

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

(23)

(24)

(25)

(26)

(27)

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

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

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

(30)

(31)

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

3.1 諧波減速器Dymola仿真模型

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

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

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

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

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

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

參 數(shù)數(shù) 值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 靜態(tài)誤差模型參數(shù)

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

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

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

圖6 諧波減速器輸入轉(zhuǎn)速Fig.6 Speed of harmonic reducer input

圖7 諧波減速器輸出轉(zhuǎn)速Fig.7 Speed of harmonic reducer output

圖8 諧波減速器動態(tài)誤差Fig.8 Dynamic error of harmonic reducer

圖9 靜態(tài)誤差和動態(tài)誤差曲線Fig.9 Static error and dynamic error curves

3.2 動態(tài)精度靈敏度及不確定性分析

Sobol敏感度分析方法[17]是一種基于方差分解的Monte Carlo方法，Sobol方法考慮了隨機輸入?yún)⒘吭谡麄€取值范圍內(nèi)對輸出響應的貢獻以及隨機參數(shù)的交互作用，Sobol敏感度指標可以直接從PCE式的系數(shù)得到。

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

(32)

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

(33)

(34)

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

(35)

其滿足

(36)

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

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

表5 動態(tài)精度多項式混沌展開式配點

表6 動態(tài)精度多項式混沌展開式系數(shù)

圖10 不同扭轉(zhuǎn)剛度k1時動態(tài)誤差變化曲線Fig.10 Dynamic error curves at different torsional stiffness k1

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

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

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

3.3 動態(tài)精度可靠性分析

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

(37)

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

項系 數(shù)項系 數(shù)項系 數(shù)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)

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

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

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

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

表8 FORM法可靠度計算結果

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

ntnfPf500140.02802000610.030540001260.031560002050.034270002540.036380003050.0381100003810.0381

4 結 論

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

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

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

3) 利用PCE方法得到動態(tài)精度的隨機統(tǒng)計特性，并將PCE方法與傳統(tǒng)的Monte Carlo方法比較，效率更高。利用PCE方法和FORM法計算一定時域內(nèi)的動態(tài)精度可靠度，比傳統(tǒng)的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，F(xiàn)AN 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] 游斌弟，趙陽.考慮非線性因素的諧波齒輪傳動動態(tài)誤差研究[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方法的翼型不確定性分析及穩(wěn)健設計[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.