劉付超，魏鵬飛，周長聰，張政，岳珠峰

西北工業(yè)大學(xué) 力學(xué)與土木建筑學(xué)院，西安 710129

運動機構(gòu)是設(shè)計用來實現(xiàn)預(yù)定運動輸出，進(jìn)而實現(xiàn)特定功能的機械裝置，在航空領(lǐng)域中有普遍的應(yīng)用。其中，平面機構(gòu)各構(gòu)件之間必須通過連接以實現(xiàn)運動傳遞。不可避免的是，伴隨著桿件連接處的不斷磨損或腐蝕，構(gòu)件連接間隙不斷擴大，此外，桿件在加工過程中也會出現(xiàn)尺寸誤差，種種因素導(dǎo)致機構(gòu)實際運動輸出與理想運動輸出并不完全切合，即會產(chǎn)生一定程度的運動誤差。研究運動誤差是否始終不超過誤差極限的問題即稱為“動態(tài)可靠性分析”[1]。其中包含兩大方面，一是在區(qū)間內(nèi)各點或整個區(qū)間段上對運動機構(gòu)的可靠性進(jìn)行分析；其次是確定影響機構(gòu)運動失效的誤差來源，即靈敏度分析。對于區(qū)間內(nèi)各點的可靠性分析，目前已經(jīng)發(fā)展出許多成熟的算法，例如一階可靠性方法[2-3]、二階可靠性方法[4]、重要性抽樣方法[5]、線性抽樣方法[6]，都可以解決類似問題。對于整個區(qū)間的時變可靠性分析，Rice首先提出了首超法[7],并取得了廣泛的應(yīng)用[8-11];其次，杜小平通過對誤差函數(shù)的一階近似估計和包絡(luò)函數(shù)法，發(fā)展出一種數(shù)值近似方法[12],使得計算精度大大提高。其他時變可靠性分析方法還包括密度函數(shù)演變法[13]、最大熵法[14]，等等。

可靠性靈敏度分析包括局部靈敏度分析和全局靈敏度分析。局部靈敏度分析為失效概率對隨機輸入變量分布參數(shù)的偏導(dǎo)數(shù)，意為隨機輸入變量分布參數(shù)的波動對于失效概率的影響程度。對于靜態(tài)結(jié)構(gòu)，文獻(xiàn)[15-17]已有成熟的局部靈敏度計算方法；對于運動機構(gòu)，一點處的局部靈敏度分析已經(jīng)實現(xiàn)[18]，更進(jìn)一步，魏鵬飛等完成了全區(qū)間段內(nèi)局部靈敏度分析[19]。與之相對應(yīng)，全局靈敏度分析旨在確定任意隨機輸入變量對于輸出結(jié)果不確定性的影響[20-23]。經(jīng)過不斷的發(fā)展，許多成熟的全局靈敏度分析方法得到了實現(xiàn)，例如基于方差的全局靈敏度分析方法[24]，矩獨立全局靈敏度分析方法[25-26]，以及基于抽樣的全局靈敏度分析方法[27]等等。

就目前所知，當(dāng)前還沒有對含旋轉(zhuǎn)鉸間隙平面運動機構(gòu)進(jìn)行全局靈敏度分析的研究工作。然而旋轉(zhuǎn)鉸間隙幾乎存在于所有運動機構(gòu)，并且對機構(gòu)的運動精度產(chǎn)生不可忽視的影響。對含旋轉(zhuǎn)鉸間隙平面運動機構(gòu)進(jìn)行全局靈敏度分析，可以有效地提高機構(gòu)可靠性，具有重要的研究意義和應(yīng)用價值。本文首先明確該種運動機構(gòu)的全局靈敏度指標(biāo)定義，其次結(jié)合相關(guān)運動機構(gòu)模型，構(gòu)造其運動誤差的包絡(luò)函數(shù)，在此基礎(chǔ)上建立含旋轉(zhuǎn)鉸間隙平面運動機構(gòu)的全局靈敏度計算方法。最后，通過兩個具體算例驗證所提方法的有效性和可行性。

1 全局靈敏度指標(biāo)定義

假定運動機構(gòu)誤差函數(shù)為g(X,θ)，X=(X1,X2,…,Xn)為機構(gòu)參數(shù)輸入變量，θ為時間因素輸入變量,誤差極限為ε。則運動機構(gòu)安全域定義為

S=X:g(X,θ)≤ε

(1)

定義安全域S指示函數(shù)為

(2)

對于指示函數(shù)IS的方差可分解為[21]

∑i≠j≠kVijk+…+V1,2,…,n

(3)

從中分離出主方差貢獻(xiàn)Vi，二階方差貢獻(xiàn)Vij，以及總方差貢獻(xiàn)VTi。

主方差貢獻(xiàn)定義為

Vi=V(IS)-EV(ISXi)

(4)

式中：EV(ISXi)為變量Xi固定后，其他變量對于指示函數(shù)方差的貢獻(xiàn)，則Vi可理解為變量Xi固定后，方差V(IS)的減小量。

二階方差貢獻(xiàn)定義為

Vij=VE(ISXi,Xj)-Vi-Vj

(5)

式中：VE(ISXi,Xj)為變量Xi、Xj對方差V(IS)的整體貢獻(xiàn)，則Vij為變量Xi、Xj兩者之間的相互作用對整體方差V(IS)的貢獻(xiàn)。

總方差貢獻(xiàn)定義為

V1,2,…,n=EV(ISX:i)

(6)

式中：X:i為變量X中除去Xi以外的n-1維向量，進(jìn)而，VTi可理解為除變量Xi外，其他變量全部固定后，變量Xi對整體方差V(IS)的貢獻(xiàn)。

進(jìn)而，通過以上3種方差貢獻(xiàn)定義出相應(yīng)的全局靈敏度指標(biāo)：

(7)

(8)

(9)

為了考察高階方差貢獻(xiàn)，利用以上3種指標(biāo)對高階指標(biāo)SHi進(jìn)行定義：

(10)

2 全局靈敏度方法

2.1 含旋轉(zhuǎn)鉸間隙平面運動機構(gòu)模型及運動誤差函數(shù)

Y=Y(X,θ)

(11)

圖1 連接間隙Fig.1 Joint clearance

以Y0表示理想運動輸出，則運動誤差函數(shù)可表示為

g(X,θ)=Y(X,θ)-Y0(θ)

(12)

2.2 基于包絡(luò)函數(shù)的可靠性分析

機構(gòu)區(qū)間可靠性定義為：對于任意θ屬于θs至θe，運動誤差函數(shù)g(X,θ)的絕對值不大于誤差極限ε的概率，即

R(θs,θe)=Pr{max(g(X,θ))≤ε∩

ming(X,θ)≥-ε}

(13)

對誤差函數(shù)g(X,θ)進(jìn)行進(jìn)一步簡化[28]：

g(X,θ)≈H(U,θ)=b0(μX,θ)+b(μX,θ)·U

(14)

式中：

(15)

并且

(16)

μgX,θ≈b0μX,θ

(17)

(18)

定義指示函數(shù)為

(19)

進(jìn)而可靠性表達(dá)式可簡化為

R(θs,θe)=Prmaxs(θ)H(U,θ)≤ε

(20)

考察機構(gòu)在某一區(qū)間內(nèi)的可靠性，可選取區(qū)間內(nèi)某些最易失效點組成聯(lián)合密度函數(shù)，進(jìn)而對其進(jìn)行積分得出機構(gòu)可靠性指標(biāo)，此即為包絡(luò)函數(shù)法。

構(gòu)造包絡(luò)函數(shù)G(X),使得

(21)

求解包絡(luò)函數(shù)G(X)，得出最可能失效點θi(i=1,2,…,t)，進(jìn)而構(gòu)造出誤差函數(shù)的均值向量為

μH=(μg1,μg2,…,μgt)

(22)

以及協(xié)方差矩陣：

(23)

其中，μgii=1,2,…,t為誤差函數(shù)在各失效點θi處的均值。Covij為g(X,θi)和g(X,θj)的協(xié)方差，i,j=1,2,…,t，具體形式為

(24)

含旋轉(zhuǎn)鉸間隙平面運動機構(gòu)可靠性表達(dá)式為

(25)

2.3 含旋轉(zhuǎn)鉸間隙平面運動機構(gòu)全局靈敏度分析

在第1節(jié)中已經(jīng)對主指標(biāo)、高階指標(biāo)以及總指標(biāo)進(jìn)行了定義。因為X為隨機向量，則IS(X)為隨機變量，且由式(2)可知，IS(X)服從兩點分布B1,R，即IS(X)等于1的概率為R，等于0的概率為1-R，因此，IS(X)的方差可以表示為

V(IS)=R(1-R)

(26)

為進(jìn)一步方便對各全局靈敏度指標(biāo)的計算，對Vi、Vij、VTi做推導(dǎo)[19]：

Vi=VE(IS|Xi)=EE2(IS|Xi)-R2

(27)

Vij=VE(IS|Xi,Xj)-Vi-Vj=

EE2(IS|Xi,Xj)-R2-Vi-Vj

(28)

VTi=EV(IS|X:i)=V(IS)-

V[E(IS|X:i)]=R-EE2(IS|X:i)

(29)

故只需求得EE2(ISXi),E[E2(ISXi,Xj)],E[E2(ISX:i)]，便可得到相應(yīng)的全局靈敏度指標(biāo)。

(30)

式中：μHMi和ΣHMi分別為HMi的均值向量和協(xié)方差矩陣。

(31)

式中：μHij和ΣHij分別為Hij的均值向量和協(xié)方差矩陣。

(32)

式中：μHTi和ΣHTi分別為HTi的均值向量和協(xié)方差矩陣。

進(jìn)而，將式(27)～式(29)和式(30)～式(32)相結(jié)合，代入式(7)～式(9)中，可分別求得全局靈敏度主指標(biāo)、高階指標(biāo)和總指標(biāo)。公式的具體推導(dǎo)過程可參考文獻(xiàn)[19]。

3 算 例

3.1 四連桿機構(gòu)

四連桿機構(gòu)在航空領(lǐng)域有極大的應(yīng)用[29]，考慮某四連桿機構(gòu)如圖2所示，1桿為驅(qū)動桿，1桿、3桿與水平線夾角θ、ψ分別為輸入變量、輸出變量。桿長服從正態(tài)分布，均值μ=(52.2,104.9,67.6,100)，標(biāo)準(zhǔn)差均為0.1；間隙尺寸(xj,yj)，j=1,2,3,4，在半徑為0.02的圓內(nèi)服從均勻分布，理想輸出函數(shù)

(33)

式中：f(x)=sinx，x∈[0, 90°]；初始輸入、輸出變量分別為θ0=95.1°、ψ0=90.6°，變化范圍分別為Δθ=120°、Δψ=60°。誤差極限ε取0.23°和0.25°?？疾?根桿L1、L2、L3、L4以及4個連接間隙C1、C2、C3、C4的全局靈敏度。以X1,X2,…,X8表示為X1=L1,X2=L2,X3=L3,X4=L4,X5=C1,X6=C2,X7=C3,X8=C4。

根據(jù)圖2可建立等式

(34)

消去ψ可得機構(gòu)實際運動輸出方程為

(35)

當(dāng)誤差極限取0.23°和0.25°時，主指標(biāo)、總指標(biāo)、高階指標(biāo)計算結(jié)果分別展示于圖3和圖4當(dāng)中，并與蒙特卡羅法計算結(jié)果相比較，蒙特卡羅法樣本點數(shù)為106。

當(dāng)誤差極限取0.23°時，由蒙特卡羅法計算出的失效概率為：Pf=0.076 004。蒙特卡羅法程序運行時間為382 s，而本文所提方法程序運行時間為64 s，原誤差函數(shù)調(diào)用次數(shù)為245次，各項指標(biāo)計算結(jié)果如圖3所示，其中，S,ST,SH分別代表主指標(biāo)、總指標(biāo)和高階指標(biāo)，S1,S2,…,S8分別為各輸入變量主指標(biāo)，ST1,ST2,…,ST8分別為各輸入變量總指標(biāo)，SH1,SH2,…,SH8分別為各輸入變量高階指標(biāo)。由圖3(a)可知，在所考察的8個變量中，X2和X4的主指標(biāo)明顯高于其他各因素,由此說明，通過降低X2與X4的不確定性，可以使得結(jié)構(gòu)的可靠性顯著提升。相反，其他各因素主指標(biāo)對結(jié)構(gòu)可靠性方差貢獻(xiàn)普遍微弱。由圖3(b)可以看出，相比于主指標(biāo)計算結(jié)果，各因素總指標(biāo)均有大幅度提升，說明各因素間相互作用明顯存在，其中，尤以X2和X4顯著。通過圖3(b)與圖3(c)對比可以發(fā)現(xiàn)，X5、X6、X7、X8的高階指標(biāo)與總指標(biāo)存在著微弱的差距，這也從另一個方面表明X5、X6、X7、X8的主指標(biāo)與二階指標(biāo)幾乎為零。

圖2 四連桿機構(gòu)Fig.2 Four-bar linkage

當(dāng)誤差極限取為0.25°時，由蒙特卡羅法計算出的失效概率為：Pf=0.0 256 045，各項指標(biāo)計算結(jié)果如圖4所示。由圖3和圖4可以看出，隨著誤差極限的增加、失效概率的下降，主指標(biāo)普遍下降、總指標(biāo)和高階指標(biāo)均微弱增長。由此說明，隨著機構(gòu)可靠性的增加，每個因素對于時變可靠性的單獨影響減小，整體影響以及各因素間相互影響普遍增加。通過對比圖3和圖4同樣可以發(fā)現(xiàn)，機構(gòu)可靠性的變化對各輸入變量的靈敏度指標(biāo)排序并沒有影響。計算成本方面，蒙特卡羅法程序運行時間為382 s，而信封函數(shù)法程序運行時間為66 s，且原誤差函數(shù)調(diào)用次數(shù)為245次，表明信封函數(shù)法計算效率高，成本小。

圖3 誤差極限為0.23°時四連桿機構(gòu)全局靈敏度分析結(jié)果Fig.3 Results of global sensitivity analysis of four-bar linkage when error limit is 0.23°

圖4 誤差極限為0.25°時四連桿機構(gòu)全局靈敏度分析結(jié)果Fig.4 Analysis of global sensitivity analysis of four-bar linkage when error limit is 0.25°

3.2 某轉(zhuǎn)向機構(gòu)

考慮某轉(zhuǎn)向機構(gòu)如圖5所示，Lb代表左側(cè)1操縱臂和右側(cè)4操縱臂的長度，Ls代表左側(cè)2轉(zhuǎn)向連接桿和右側(cè)3轉(zhuǎn)向連接桿的長度，W表示下側(cè)輪軸長度，H為輪軸與齒條軸5之間的距離。以齒條軸5的橫向位移D為輸入變量，左側(cè)操縱臂與水平線夾角θ1為輸出變量。X=(X1,X2,X3,X4)=(Lb,Ls,W,H)代表4個隨機尺寸變量，均值μ=(111,283.5,650.24,83.5) 10-3m，標(biāo)準(zhǔn)差為σ=(0.1,0.1,0.1,0.1) 10-3m。間隙變量(xj,yj),j=1,2，在半徑為0.1的圓內(nèi)服從均勻分布。機構(gòu)理想運動輸出假定為實際運動輸出函數(shù)在尺寸變量和間隙變量均值處的取值?？疾?個尺寸因素Lb、Ls、W、H以及2個間隙因素C1、C2的全局靈敏度。分別以X1,X2,…,X6表示為:X1=Lb,X2=Ls,X3=W,X4=H,X5=C1,X6=C2。

根據(jù)圖5可得等式：

(36)

消去θ2可得機構(gòu)實際運動輸出方程為

(37)

當(dāng)誤差極限取0.4°時，主指標(biāo)、總指標(biāo)、高階指標(biāo)計算結(jié)果展示于圖6當(dāng)中。蒙特卡羅法樣本點數(shù)為106。

當(dāng)誤差極限取0.4°時，由蒙特卡羅法計算出的失效概率為：Pf=0.0 663 576。蒙特卡羅法程序運行時間為248 s，信封函數(shù)法程序運行時間為43 s，原誤差函數(shù)調(diào)用次數(shù)為245次。由圖6(a)可知，所考察的6個變量的主指標(biāo)普遍較小，由此也導(dǎo)致兩種方法的計算結(jié)果在圖像上存在明顯的相對差異，但通過比較可知，絕對偏差較小，所提方法依舊具有足夠的準(zhǔn)確度。從圖中可以看出，相比于其他4個因素，X1和X2具有較高的主指標(biāo)，由此說明，對各輸入變量均降低一定程度的不確定性，相對于其他輸入變量，X1和X2可使結(jié)構(gòu)的可靠度得到更多的提高。由圖6(b)可以看出，各因素總指標(biāo)相對于主指標(biāo)均有大幅度提升，這也使得兩種方法的計算結(jié)果具有更高的貼合度，同時說明各輸入變量間具有明顯的相互作用。在全部輸入變量中，X1、X2具有較高的總指標(biāo)，其次為X5、X6，值得注意的是，兩個間隙因素X5、X6具有相同的總指標(biāo)，說明兩因素整體上對機構(gòu)可靠性具有相同影響。在圖6(b)與圖6(c)的對比中可以發(fā)現(xiàn)，各輸入變量高階指標(biāo)與總指標(biāo)間具有明顯的差距，且由圖6(a)可知各變量主指標(biāo)較小，從而表明各輸入變量兩兩之間具有顯著的相互作用。

圖5 某轉(zhuǎn)向機構(gòu)示意圖Fig.5 Diagram of a steering gear

圖6 某轉(zhuǎn)向機構(gòu)全局靈敏度分析結(jié)果Fig.6 Results of global sensitivity analysis of steering gear

4 結(jié) 論

1) 對于含旋轉(zhuǎn)鉸間隙的平面運動機構(gòu)，本文在時變運動機構(gòu)可靠性靈敏度分析方法的基礎(chǔ)上，結(jié)合包絡(luò)函數(shù)法，發(fā)展出一種高效的計算含旋轉(zhuǎn)鉸間隙平面運動機構(gòu)全局靈敏度的分析方法。

2) 將本文所發(fā)展的高效計算方法應(yīng)用到兩個實際算例中，并與蒙特卡羅法相對比，結(jié)果顯示，本文所提方法所得結(jié)果與蒙特卡羅法所得結(jié)果非常接近，程序運行時間大大縮短，且誤差函數(shù)調(diào)用次數(shù)較少，表明本文所提方法精確度高，計算成本小，具有良好的工程應(yīng)用價值。