柔性機構時變可靠性分析的時變隨機響應面法

闞琳潔，張建國，邱繼偉

(1. 北京航空航天大學 可靠性與系統工程學院, 北京 100191；2. 北京航空航天大學 可靠性與環境工程技術重點實驗室, 北京 100191)

在航空航天、機器人、汽車、艦船等諸多行業中，柔性機構都具有舉足輕重的地位，近幾年來備受人們的關注[1-2]。在大多數實際工程問題中，柔性機構的系統參數如材料性能、幾何形狀、剛度、阻尼系數、制造公差等都具有隨機性，使得柔性機構可靠性問題尤為重要。對于柔性機構時變可靠性分析問題，傳統且有效的方法是利用經典的結構時變可靠性理論建立極限狀態函數，采用Monte Carlo仿真分析(MCS) 、解析法[3]或上穿率法[4]。然而，由于其系統動力學方程的非線性、耦合、時變特征，一次模型的計算量通常是相當高的，采用MCS法計算效率會很低。而且極限狀態函數在非線性影響下難以顯示化表達，使得解析法和上穿率法也難以實施[5]。因此，一些學者開始研究基于代理模型的柔性機構時變可靠性分析技術，并取得了一定的成果，例如響應面[6]、神經網絡[7]、支持向量機[8]、混沌多項式(Polynomial Chaos Expansion，PCE)[9]等。

由Wiener[10]提出的混沌多項式是一種正交展開方法，標準正交函數系由Hermit基函數構造。Liu等[11]對其進行了擴展，提出了基于Askey-scheme正交基的廣義混沌多項式(Generalized Polynomial Chaos，GPC)，常用于隨機系統的建模、不確定性傳播和量化、以及求解任意概率分布的隨機微分方程等方面。作為新興的代理模型技術，GPC屬于隨機數學模型。相較于其他代理模型，GPC不受擬合點位置的影響，近似精度高，并且可以獲得系統響應精確解的全部信息[12]，但應用于柔性機構時變可靠性分析的研究成果較少。Guo等[13]利用GPC的非侵入式算法結合模糊概率理論，提出了考慮固有-認知混合不確定信息的機構運動可靠度分析方法。并利用移動最小二乘法來近似表達時變GPC系數的變化趨勢，提出了考慮磨損的機構時變全局靈敏度分析法，以及考慮磨損的機構運動的時變可靠性優化設計方法[14]。趙寬[15]利用基于PCE的隨機響應面法，對隨機參數旋轉柔性梁和雙連桿機械臂進行了機構性能(運動、強度、剛度)的靜態可靠性分析。雖然混沌多項式模型計算精度高，但是系統隨機參數數量多，剛柔耦合引起的非線性問題，以及系統結構的復雜性，會使得混沌多項式理論用于時變可靠性分析時計算成本迅速增長。

綜上所述，柔性機構的系統動力學模型一般為具有時變、高非線性、強耦合特征的微分方程組，其系統可靠性問題同時具有時變耦合的特點。傳統的柔性機構可靠性分析方法不僅難以描述柔性機構的時變耦合特征，而且模型精度差、計算效率低。為了提高計算效率，描述柔性機構在剛柔耦合影響下的時變可靠性特征，本文提出基于模態綜合法(Component Mode Synthesis，CMS)和GPC的時變隨機響應面法。建立了柔性機構時變可靠性分析的隨機代理模型，并用MCS法進行系統運動時變可靠度求解。

1 柔性機構動力學的隨機場離散和模型降階

1.1 隨機場離散

考慮機構系統的材料特性與幾何尺寸的隨機性，基于混合坐標法，利用有限元法對柔性體的變形進行描述，結合拉格朗日方程可推導出如下形式的柔性機構系統動力學的隨機模型[16]

(1)

如前所述，柔性機構的材料特性參數具有不確定性，并在空間域中連續變化，可用隨機場進行描述。設具有協方差函數C(x1,x2)的隨機場可表示為h(x,θ)，用Karhunen-Loeve展開法進行隨機場離散可得[17]

(2)

(3)

則有

(4)

1.2 動力學模型降階的CMS法

采用有限元法對柔性體變形進行描述時，坐標數目龐大，動力學模型式(1)為非線性、耦合的微分方程組。為了提高計算效率，本文采用CMS法對剛柔耦合動力學模型進行降階，并且假設柔性體模態為確定性模態[18]。在模型降階過程中暫不考慮剛體運動的影響，僅對柔性體的彈性變形進行降階，其動力學方程為

(5)

(6)

(7)

(8)

(9)

降階后的柔性體運動方程可表示為

(10)

(11)

(12)

(13)

式中:Ξ∈RD×H為正交化的CMS模態矩陣。利用均值處的模態矩陣對隨機動力學模型式(4)進行降階，則廣義坐標q的正交投影變換可以寫為

(14)

(15)

2 基于CMS-GPC的時變隨機響應面

2.1 時變隨機響應面

利用GPC模型可以有效地對系統響應的隨機特征進行描述，并具有顯示化和線性化的作用，但是GPC對系統模型也產生了一定的擴階影響，廣義坐標q的自由度從RC擴展到RC·(M+1)。為了提高計算效率，本文將CMS方法與GPC方法相結合，提出時變隨機響應面的隨機代理模型。

(16)

式中:t∈[0,T];{Ψj(ξ),j=0,1,2,,M}為隨機函數空間中的Askey-scheme正交基，且Ψ∈R(6+H)×[(6+H)·(M+1)]，ξ=(ξi1,ξi2,,ξin)為N維隨機向量，a(t)∈R(6+H)·(M+1)為多項式時變系數，M為截斷系數。若Askey-scheme正交基的最高階次為p，則最高階次p、隨機向量維數N和截斷系數M滿足如下等式關系

(17)

對應不同分布類型的Askey-scheme正交基如表1所示。

表1 對應不同分布類型的Askey-scheme正交基Tab.1 Askey-scheme for selecting polynomials corresponding to certain types of distributions

任意兩個Askey-scheme正交基之間滿足

(18)

(19)

其次，根據式(14)和式(16)，系統廣義坐標q可以寫為兩級正交展開變換形式

(20)

也可以寫為如下截斷形式

(21)

(22)

原GPC模型計算問題的大小從C·(M+1)下降到(6+H)·(M+1)，其中(6+H)<

2.2 時變隨機響應面的時變系數求解

用Λ(ξ)對式(4)進行Galerkin投影，可得

整理可得

(23)

則式(23)為M+1個確定性的二階微分方程，可求得時變隨機響應面的時變系數a(t)。

(24)

3 柔性機構時變可靠性分析

3.1 時變可靠度計算方法

(25)

(26)

在時間[0,T]內，運動精度的時變累積失效概率為

Pf,c(0,T)=Pr(?t∈[0,T]|G(ξ,t)≤0)

(27)

在建立柔性機構時變可靠性模型的基礎上，采用MCS法進行時變可靠性分析。將時間離散成N個時間點，每個時間間隔為Δt=T/(N-1)，對應的響應面G(ξ,t)離散成{G(ξ,ti),ti=(i-1)Δt且i=1,2,,N}，則在時間[0,ti]內運動精度的時變累積失效概率為[19]

(28)

式中:NMCS為總的仿真次數;Nj,fature為在時間[ti,ti+1]內系統失效的次數。

3.2 柔性機構時變可靠性分析流程

歸納上述基于CMS-GPC的柔性機構時變可靠性建模分析過程，給出如下分析步驟:

步驟1確定系統隨機參數的分布類型。對高斯隨機場可采用可以通過 Karhunen-Loeve 展開轉化成離散的標準隨機向量。

步驟2建立柔性機構動力學模型，依據步驟1對隨機參數的處理，將動力學模型進行隨機展開，形式如式(2)，并確定模型的初始條件。

步驟3求解模態降階矩陣Φ。采用CMS法對剛柔耦合動力學模型進行降階，并且假設柔性體模態為確定性模態，得到均值處的模態降階矩陣Φ。

步驟4確定GPC模型中正交基Ψ的階次p和截斷系數M，建立系統響應的時變隨機響應面q(ξ,t)=ΦΨ(ξ)a(t)=Λ(ξ)a(t)。利用Galerkin投影法求解響應面的時變系數，分析系統響應的時變均值和時變方差。

步驟5基于時變隨機響應面，建立系統運動精度的時變可靠性模型，并利用MCS法求解時變可靠度。

4 案例分析

以兩連桿柔性機械臂為例，對本文提出的方法進行驗證。圖1為兩連桿柔性機械臂示意圖，在空間零重力環境下做平面大范圍運動，初始狀態為水平位置，無變形，無初始速度和加速度。

圖1 兩連桿柔性機械臂Fig.1 Two-link flexible manipulator

組成機械臂的桿1和桿2為均質的Euler-Bernoulli梁，橫截面均為矩形，服從正態分布，相關參數由表2給出。臂桿長度為L1=L2=1.5 m，作用在臂桿上的驅動力矩分別為τ1(t)=60sin2(2πt)+25，τ2(t)=50sin2(2πt+0.5)+15。關節處集中質量為m1=10 kg，桿2末端負載集中質量為m2=5 kg。機械臂的材料參數服從正態分布，相關參數由表3給出。

表2 兩連桿機械臂橫截面參數Tab.2 Section sizes of two-link flexible arm

表3 兩連桿機械臂材料參數Tab.3 Material parameters of two-link flexible arm

兩臂桿各離散5個單元，建立兩連桿柔性機械臂動力學模型，并利用式(2)～式(4)進行隨機場離散。其次，桿1和桿2分別取前兩階模態，建立兩連桿柔性機械臂末端X方向位移x(ξ,t)和Y方向位移y(ξ,t)的時變隨機響應面模型，取Λ(ξ)的最高階次為4階，根據式(17)則響應面的截斷系數為210。

(29)

圖2和圖3分別給出了兩連桿柔性機械臂末端X方向位移的均值和均方差隨時間變化的規律，并與Monte-Carlo模擬105次的計算結果進行對比，時變均值的最大相對誤差為0.17%，時變方差的最大相對誤差為0.56%，誤差較小，具有較高的擬合精度。

圖2 機械臂末端X方向位移均值Fig.2 Mean displacement of the end node in X direction

圖3 機械臂末端X方向位移方差Fig.3 Displacement covariance of the end node in X direction

圖4和圖5分別給出了兩連桿柔性機械臂末端Y方向位移y(ξ,t)的均值和均方差隨時間變化的規律。同樣與Monte-Carlo模擬105次的結果相比，時變均值的最大相對誤差為0.17%，時變方差的最大相對誤差為0.56%，同樣誤差較小，說明本文所提出的時變耦合響應面具有較高的計算精度。同時，計算規模與僅用GPC模型相比，減小了80%。

圖4 機械臂末端Y方向位移均值Fig.4 Mean displacement of the end node in Y direction

(30)

Y方向運動精度的時變極限狀態函數為

(31)

在此基礎上，采用MCS法計算計算X方向和Y方向定位精度的時變累積失效概率。圖6和圖8分別為X方向和Y方向的累積失效概率在0～2 s內的變化趨勢，并且與MCS仿真105次進行了結果對比，其相對誤差隨時間的變化趨勢如圖7和圖9所示。從結果來看，與僅用MCS法相比，累積失效概率的變化趨勢相同。其中X方向累積失效概率的最大相對誤差為0.75%，Y方向累積失效概率的最大相對誤差為0.76%，故誤差結果在可接受范圍內，說明本文提出的方法是可行的。

圖6 機械臂末端X方向定位精度累積失效概率圖7 末端X方向累積失效概率的相對誤差圖8 機械臂末端Y方向定位精度累積失效概率圖9 末端Y方向累積失效概率的相對誤差Fig.6 The cumulative failure probability of the positioning accuracy of the end node in X directionFig.7 The relative error of the cumulative failure probability of the end node in X directionFig.8 The cumulative failure probability of the positioning accuracy of the end node in Y directionFig.9 The relative error of the cumulative failure probability of the end node in Y direction

5 結 論

(1)本文針對柔性機構系統時變可靠性分析的“時變耦合”問題，提出了基于CMS-GPC的時變隨機響應面模型，在保證擬合精度的同時，降低了計算問題的規模，提高了效率。其次，與傳統時變可靠性分析的響應面模型不同的是，本文所建立的時變隨機響應面模型是描述系統隨機響應特征的，而非對特定時間下極限狀態函數的擬合。

(2)在時變隨機響應面基礎上，可直接建立顯示化的柔性機構運動時變可靠性模型，并給出了計算運動精度時變累積失效概率的MCS方法。與傳統時變可靠性分析方法相比，本方法不需要在設計驗算點處對特定時間下的極限狀態函數重復建立響應面，因此提高了計算效率。

(3)通過兩連桿柔性機械臂的案例對本文的方法進行了可行性和有效性驗證，計算精度與MCS法相差無幾，對工程應用具有一定的指導意義。