劉延浩 倪冰雨 田萬一 姜潮

湖南大學(xué)機(jī)械與運(yùn)載工程學(xué)院，長沙，410006

摘要：

針對動態(tài)不確定性激勵作用下的剛?cè)狁詈舷到y(tǒng)動力學(xué)問題，提出了一種基于區(qū)間過程模型的動力學(xué)不確定性序列模擬方法，旨在通過區(qū)間過程序列抽樣及剛?cè)狁詈蟿恿W(xué)序列模擬計算結(jié)構(gòu)振動與機(jī)構(gòu)運(yùn)動等系統(tǒng)動態(tài)響應(yīng)的上下邊界。介紹了中心剛體柔性梁剛?cè)狁詈舷到y(tǒng)動力學(xué)方程的構(gòu)建與數(shù)值求解方法。針對動態(tài)不確定性激勵作用下的剛?cè)狁詈舷到y(tǒng)動力學(xué)分析，引入?yún)^(qū)間過程模型及其區(qū)間K-L展開對動態(tài)不確定性進(jìn)行了度量和高效表征，提出了一種求解系統(tǒng)機(jī)構(gòu)運(yùn)動與結(jié)構(gòu)振動等動力學(xué)響應(yīng)上下邊界的序列模擬方法。該方法利用序列模擬策略識別當(dāng)前模擬序列中對動態(tài)響應(yīng)上邊界或下邊界具有貢獻(xiàn)的區(qū)間過程參數(shù)樣本集，作為下一模擬序列中的局部加密抽樣中心，可有效避免直接蒙特卡羅模擬在計算動力學(xué)響應(yīng)上下邊界時因過多無效抽樣模擬而導(dǎo)致的低效收斂問題。最后，通過三個算例對所提方法的有效性進(jìn)行了驗(yàn)證。研究結(jié)果表明，對于剛?cè)狁詈舷到y(tǒng)大范圍運(yùn)動及振動響應(yīng)上下邊界的求解，序列模擬方法相比于直接蒙特卡羅模擬方法具有更好的計算效率和計算精度。

關(guān)鍵詞：動態(tài)不確定性；區(qū)間過程模型；區(qū)間K-L展開；剛?cè)狁詈蟿恿W(xué)；序列模擬方法

中圖分類號：TH113

DOI：10.3969/j.issn.1004132X.2024.05.002

開放科學(xué)（資源服務(wù)）標(biāo)識碼（OSID）：

A Sequential Simulation Method for Dynamic Uncertainty Analysis of

Rigid-flexible Coupling Systems under Interval Process Excitations

LIU Yanhao? NI Bingyu? TIAN Wanyi? JIANG Chao

School of Mechanical and Vehicle Engineering，Hunan University，Changsha，410006

Abstract： For the dynamic problem of rigid-flexible coupling systems under dynamic uncertain excitations， an interval process model-based sequential simulation method was proposed for uncertainty analysis， which aimed to obtain the upper and lower bounds of the system dynamic responses such as structural vibrations and mechanism kinematics， by sequential sampling of the interval process and the rigid-flexible coupling dynamics simulations. The construction and numerical solution of the dynamic equation of the rigid-flexible coupling systems with central rigid body and flexible beam were introduced. Aiming at the dynamic analysis of rigid-flexible coupling systems under uncertain dynamic excitations， the interval process model and the interval K-L expansion were introduced to quantify and represent the dynamic uncertainty efficiently， and a sequential simulation method was proposed to solve the upper and lower bounds of the dynamic responses of the system mechanism motions and structural vibrations. The method used a sequential simulation strategy to identify the interval process parameter sample sets that contributed to the upper or lower bounds of dynamic responses in the cur rent simulation sequence， and served as the local encrypted sampling center in the next simulation sequence， which might effectively avoid the inefficient convergence problem caused by excessive invalid sampling simulations when calculating the upper and lower bounds of dynamic response in direct Monte Carlo simulation. Finally， three examples were given to verify the effectiveness of the proposed method. The results show that the sequential simulation method has better computational efficiency and accuracy than that of the direct Monte Carlo simulation method for solving the upper and lower bounds of the rigid-flexible coupling systems large overall motions and vibration responses.

Key words： dynamic uncertainty； interval process model； interval K-L expansion； rigid-flexible coupling dynamics； sequential simulation method

收稿日期：20240415

基金項(xiàng)目：國家自然科學(xué)基金（52175224，52235005，52075156）

0? 引言

柔性多體系統(tǒng)動力學(xué)研究由剛體和柔性體組成的復(fù)雜機(jī)械系統(tǒng)在經(jīng)歷大范圍空間運(yùn)動時的動力學(xué)行為，在航天器、機(jī)器人、兵器與車輛等工程領(lǐng)域得到了廣泛關(guān)注。剛?cè)狁詈舷到y(tǒng)的動力學(xué)問題是柔性多體系統(tǒng)動力學(xué)的本質(zhì)，主要研究柔性體變形與其大范圍空間運(yùn)動之間的相互耦合作用，以及這種耦合所導(dǎo)致的動力學(xué)效應(yīng)［1-3］。然而，由于服役環(huán)境的復(fù)雜性和運(yùn)行工況的多樣性，系統(tǒng)通常受到不確定性因素影響導(dǎo)致機(jī)構(gòu)運(yùn)動、振動等動力學(xué)響應(yīng)產(chǎn)生波動。例如，高速與高精度裝配機(jī)器人的定位精度必須考慮機(jī)械臂在高速運(yùn)動中與其自身彈性變形之間的剛?cè)狁詈闲?yīng)，然而基體振動、外部擾動等動態(tài)不確定性的影響將降低機(jī)器人的裝配精度。因此，合理量化上述動態(tài)不確定性參數(shù)對剛?cè)狁詈蟿恿ο到y(tǒng)的影響，精確預(yù)測系統(tǒng)動力學(xué)響應(yīng)的波動特性，并在設(shè)計中加以充分考慮，對提高系統(tǒng)性能及可靠性具有重要意義。

區(qū)間過程激勵下剛?cè)狁詈舷到y(tǒng)動態(tài)不確定性分析的序列模擬方法——劉延浩? 倪冰雨? 田萬一等

中國機(jī)械工程 第35卷 第5期 2024年5月

目前，在剛?cè)狁詈舷到y(tǒng)動力學(xué)不確定性分析領(lǐng)域，國內(nèi)外學(xué)者已經(jīng)進(jìn)行了相關(guān)研究。吳勝寶等［4］對自由大范圍運(yùn)動情況下剛體柔性梁系統(tǒng)的剛?cè)狁詈蟿恿W(xué)問題進(jìn)行了研究，分析了系統(tǒng)在不同外界激勵作用下的耦合動力學(xué)響應(yīng)。為了減少工業(yè)機(jī)器人的制造和各種環(huán)境等隨機(jī)不確定性因素造成結(jié)構(gòu)系統(tǒng)性能下降，胡啟國等［5］采用一次二階矩法得到臂部可靠性指標(biāo)，并用可靠性指標(biāo)反映不確定性因素的影響，并提出了一種考慮剛?cè)狁詈系墓I(yè)機(jī)器人多目標(biāo)可靠性拓?fù)鋬?yōu)化方法。辛鵬飛等［6］以柔性空間機(jī)械臂為例，運(yùn)用基于Chebyshev多項(xiàng)式的區(qū)間擴(kuò)張函數(shù)，將含區(qū)間參數(shù)的微分代數(shù)方程轉(zhuǎn)化為Chebyshev多項(xiàng)式插值點(diǎn)處的確定參數(shù)的動力學(xué)方程，研究得到了一維區(qū)間參數(shù)和多維區(qū)間參數(shù)影響下機(jī)械臂系統(tǒng)的動力學(xué)響應(yīng)區(qū)間邊界，形成了預(yù)測機(jī)械臂末端軌跡區(qū)間的方法。WANG等［7］從不確定性和系統(tǒng)工程的角度研究了火炮系統(tǒng)結(jié)構(gòu)動力響應(yīng)的不確定優(yōu)化問題，針對火炮發(fā)射過程中內(nèi)部彈道、發(fā)射載荷等參數(shù)波動對彈丸運(yùn)動的影響，提出了一種考慮魯棒性和經(jīng)濟(jì)性的火炮結(jié)構(gòu)動力響應(yīng)區(qū)間不確定優(yōu)化方法。李勇［8］基于某火炮供藥機(jī)構(gòu)剛?cè)狁詈蟿恿W(xué)模型，對模塊藥放入、貯存和回位等多個機(jī)構(gòu)動作的不確定性進(jìn)行了仿真分析。宋清林［9］考慮運(yùn)動副間隙和協(xié)調(diào)臂柔性，研究了多個參數(shù)波動對火炮協(xié)調(diào)器動態(tài)特性的不同影響，并進(jìn)一步分析了協(xié)調(diào)運(yùn)動的可靠性。YANG等［10］將剛?cè)狁詈闲l(wèi)星中的多源不確定性量化為區(qū)間參數(shù)，基于標(biāo)稱剛?cè)狁詈闲l(wèi)星動力學(xué)模型建立了相應(yīng)的耦合狀態(tài)方程，同時，考慮代理模型中的多項(xiàng)式混沌展開，提出了一種區(qū)間分析方法預(yù)測不確定姿態(tài)振動控制系統(tǒng)的不確定狀態(tài)響應(yīng)。根據(jù)文獻(xiàn)調(diào)研，目前針對剛?cè)狁詈蟿恿W(xué)系統(tǒng)的不確定性分析主要關(guān)注加工誤差、材料屬性、條件約束等參量不確定性，尚缺乏對動態(tài)不確定性因素的考慮。然而，對于工程實(shí)際中的諸多剛?cè)狁詈蟿恿W(xué)系統(tǒng)，振動激勵、交變溫度載荷等動態(tài)不確定性對此類系統(tǒng)尤其是高速高精度執(zhí)行機(jī)構(gòu)等系統(tǒng)的動力學(xué)性能具有重要影響，相關(guān)研究對后續(xù)系統(tǒng)高性能保障和高可靠性設(shè)計具有重要的理論和應(yīng)用價值。但不同于剛體動力學(xué)系統(tǒng)和結(jié)構(gòu)振動系統(tǒng)，剛?cè)狁詈蟿恿W(xué)系統(tǒng)通常同時伴隨機(jī)構(gòu)運(yùn)動和結(jié)構(gòu)振動的復(fù)雜耦合作用，使其動力學(xué)方程呈現(xiàn)出高度非線性、強(qiáng)耦合和時變性等特點(diǎn)；在隨機(jī)動態(tài)載荷等動態(tài)不確定性作用下，對系統(tǒng)機(jī)構(gòu)運(yùn)動及結(jié)構(gòu)振動的耦合動力學(xué)響應(yīng)進(jìn)行求解屬于強(qiáng)非線性系統(tǒng)的動力學(xué)不確定性分析，這對分析方法在計算效率和精度等方面提出了挑戰(zhàn)。

本文綜合考慮剛?cè)狁詈舷到y(tǒng)動態(tài)不確定性問題的工程特性，提出了一種基于區(qū)間過程模型的動力學(xué)不確定性序列模擬方法，以獲得剛?cè)狁詈舷到y(tǒng)機(jī)構(gòu)運(yùn)動狀態(tài)參數(shù)及結(jié)構(gòu)振動響應(yīng)隨時間變化的上下邊界。該方法在剛?cè)狁詈舷到y(tǒng)動態(tài)不確定性分析方面的優(yōu)勢主要包括：①對于許多實(shí)際工程問題，由于測試條件或經(jīng)濟(jì)成本的限制等，通常外部載荷、路面激勵等動態(tài)不確定參數(shù)的測試數(shù)據(jù)量難以滿足隨機(jī)過程建模的要求，而區(qū)間過程由于在任意時刻使用區(qū)間而非精確概率分布來描述參數(shù)的變化范圍，在整個時間歷程上通過兩條邊界曲線描述時變參數(shù)的不確定性，可很大程度上降低對參數(shù)樣本量的依賴；②對于諸多柔性多體系統(tǒng)，所受到的動態(tài)不確定性（如摩擦磨損、沖擊振動等）通常來源于上游/前端間隙、磨損量、溫度等物理屬性的綜合影響，而此類前端參數(shù)的不確定性通常認(rèn)為屬于有界不確定性，因此在復(fù)雜系統(tǒng)中演化影響形成的動態(tài)不確定性也屬于有界動態(tài)不確定性，適用于描述為區(qū)間過程；③基于區(qū)間過程的剛?cè)狁詈蟿恿W(xué)不確定性求解可獲得結(jié)構(gòu)動態(tài)響應(yīng)邊界或機(jī)構(gòu)運(yùn)動軌跡包絡(luò)，這種邊界或包絡(luò)形式對機(jī)構(gòu)運(yùn)動精度分析和機(jī)構(gòu)運(yùn)動干涉判斷十分有利和方便；④剛?cè)狁詈蟿恿W(xué)問題通常具有強(qiáng)非線性，并且由于涉及大范圍機(jī)構(gòu)運(yùn)動和結(jié)構(gòu)振動耦合計算，往往計算復(fù)雜，因此，為便于工程分析，適宜采用非嵌入式不確定性求解方法?，F(xiàn)有考慮區(qū)間過程的不確定性分析主要針對結(jié)構(gòu)振動問題，且基本上都是處理線彈性振動系統(tǒng)和少數(shù)弱非線性振動系統(tǒng)。本文方法首先引入?yún)^(qū)間過程模型及區(qū)間K-L展開方法對剛?cè)狁詈蟿恿W(xué)系統(tǒng)中的動態(tài)不確定性進(jìn)行量化和高效表征，進(jìn)而利用序列模擬策略識別當(dāng)前剛?cè)狁詈蟿恿W(xué)分析模擬序列中對系統(tǒng)機(jī)構(gòu)運(yùn)動與結(jié)構(gòu)振動等動力學(xué)響應(yīng)上邊界或下邊界具有貢獻(xiàn)的區(qū)間過程參數(shù)樣本集，并作為下一模擬序列中的局部加密抽樣中心，有效避免了直接蒙特卡羅模擬在計算動力學(xué)響應(yīng)上下邊界時因過多無效抽樣模擬而導(dǎo)致的低效收斂問題，從而有效提高了不確定性求解效率。此外，本文簡要介紹了剛?cè)狁詈舷到y(tǒng)動力學(xué)建模與數(shù)值求解方法，給出了基于區(qū)間K-L展開的剛?cè)狁詈舷到y(tǒng)動力學(xué)不確定性分析序列模擬求解方法，并通過三個算例驗(yàn)證了本文方法的有效性。

1? 剛?cè)狁詈舷到y(tǒng)動力學(xué)分析簡介

1.1? 剛?cè)狁詈舷到y(tǒng)物理模型

圖1所示的中心剛體柔性梁剛?cè)狁詈舷到y(tǒng)屬于一類典型的剛?cè)狁詈蠁栴}［4］，其中，中心剛體做平面運(yùn)動，其上以懸臂方式連接柔性梁。坐標(biāo)系OXY為系統(tǒng)的慣性坐標(biāo)系，在柔性梁上建立浮動坐標(biāo)系oxy，圖1中o為剛體與柔性梁的連接處；θ為ox軸與OX軸的夾角；x為未變形時柔性梁的軸線；剛體質(zhì)量為m，剛體繞質(zhì)心c的轉(zhuǎn)動慣量為Jch，半徑為R；柔性梁的質(zhì)量為m1，長度為L，密度為ρ，橫截面積為S，彈性模量為E，截面慣性矩為I；τ為作用在轉(zhuǎn)軸上的力矩；r為柔性梁變形后一點(diǎn)z到O點(diǎn)的矢徑；ω1、ω2分別為梁的軸線縱向伸長量和為橫向彎曲變形量；ωc為耦合變形量。

1.2? 剛?cè)狁詈舷到y(tǒng)動力學(xué)方程構(gòu)建與求解

多體系統(tǒng)動力學(xué)方程的推導(dǎo)方法有多種，一般常用的有牛頓歐拉法、拉格朗日法和高斯法等［11-12］，在構(gòu)建多體系統(tǒng)動力學(xué)方程中，多采用拉格朗日法，該類方程可直接表示為系統(tǒng)控制輸入的函數(shù)。對于完整系統(tǒng)用廣義坐標(biāo)表示的動力學(xué)方程，通常系指第二類拉格朗日方程。

如圖1所示，從慣性坐標(biāo)系OXY原點(diǎn)引往柔性梁變形后一點(diǎn)z（x，y）的矢徑r為［4］

r=rc+ro+ρo+u（1）

式中，rc為剛體質(zhì)心關(guān)于慣性坐標(biāo)系的矢徑；ro為剛體質(zhì)心至浮動坐標(biāo)系原點(diǎn)的矢徑；ρo為梁變形前z點(diǎn)關(guān)于浮動坐標(biāo)系原點(diǎn)的矢徑；u為梁變形前后末端位移矢徑。

設(shè)λ為浮動坐標(biāo)系相對于慣性坐標(biāo)系的方向余弦矩陣，則矢徑r在慣性系下的坐標(biāo)陣為

r=rc+λ（ro+ρo+u）（2）

rc=（X，Y）T? ro=（R，0）T

ρo=（x，y）T? u=（ux，uy）

λ=cos θ－sin θsin θcos θ

將式（2）對時間t求導(dǎo)可得到z點(diǎn)速度為

r·=r·c+λ·（ro+ρo+u）+λu·（3）

系統(tǒng)的動能和勢能分別為

Ek=12∫Vρr·Tr·dV+12mr·2c+12Jchθ2（4）

Ep=12∫L0ES（ω1x）2dx+12∫L0EI（2ω2x2）2dx（5）

式中，V為體積。

求解單自由度柔性臂的振動模型可采用假設(shè)模態(tài)法對模態(tài)坐標(biāo)進(jìn)行解耦處理，再對結(jié)果進(jìn)行線性變換得到實(shí)際的坐標(biāo)響應(yīng)。采用假設(shè)模態(tài)法描述柔性梁的變形，縱向變形ω1和橫向變形ω2分別表示為［4］

ω1（x，t）=Φx（x）A（t）

ω2（x，t）=Φy（x）B（t）（6）

Φx（x）=（Φ（1）x（x），Φ（2）x（x），…，Φ（i）x（x），…，Φ（N）x（x））

A（t）=（A1（t），A2（t），…，Ai（t），…，AN（t））T

Φy（x）=）（Φ（1）y（x），Φ（2）y（x），…，Φ（i）y（x），…，Φ（N）y（x））

B（t）=（B1（t），B2（t），…，Bi（t），…，BN（t））T

Φ（i）x（x）=sin（2i－12πxL）? i=1，2，…，N

Φ（i）y（x）=cosh（βix）－cos（βix）－γi（sinh（βix）－

sin（βix））? i=1，2，…，N

γi=cosh（βiL）+cos（βiL）sinh（βiL）+sin（βiL）

式中，Φx（x）、Φy（x）分別為梁的縱向振動和橫向振動的模態(tài)函數(shù)行矢量;A（t）、B（t）分別為縱向振動和橫向振動的模態(tài)坐標(biāo)列矢量；N為模態(tài)函數(shù)的階數(shù)；βi為模態(tài)函數(shù)行矢量特征值，可由特征方程1+cosh（βiL）cos（βiL）=0求出。

采取廣義坐標(biāo)q=（X，Y，θ，A，B），運(yùn)用第二類拉格朗日方程可得

ddt（Tq·）－Tq=－Uq+Fq（7）

Fq=（FX，F(xiàn)Y，τ，0，0）

式中，F(xiàn)q為外主動力對應(yīng)的廣義力；FX、FY分別為剛體上合力主矢在X、Y方向上投影；τ為剛體上合力關(guān)于質(zhì)心的矢矩。

將式（4）和式（5）代入式（7），得到剛?cè)狁詈舷到y(tǒng)動力學(xué)方程為

M（q）q¨=Q（q，q·，t） （8）

其中，M、Q分別為廣義質(zhì)量陣和廣義力陣，可分別表示為

M=M110M13M14M150M22M23M24M25M31M32M33M34M35M41M42M43M440M51M52M530M55

Q=［QX? QY? Qθ? QA? QB］T（9）

進(jìn)一步地，可以得到如下方程：

M110M13M14M15

0M22M23M24M25

M31M32M33M34M35

M41M42M43M440

M51M52M530M55X¨Y¨θ¨A¨B¨=QXQYQθQAQB（10）

式（9）、式（10）中部分變量為矩陣，各變量的具體計算表達(dá)式可以參考文獻(xiàn)［4］。

目前，對上述剛?cè)狁詈舷到y(tǒng)動力學(xué)方程的求解已經(jīng)可以使用較為成熟的數(shù)值方法［13］，主要包括：四階龍格庫塔法（fourth-order Runge-Kutta method，RK4）［14］、自適應(yīng)變步長的龍格庫塔費(fèi)爾博格法（Runge-Kutta-Fehlberg method，RKF45）；三階吉爾預(yù)估校正法（third-order Gear predictor-

corrector algorithm，3 Gear）；基于紐馬克（Newmark）法［15］的希爾伯特修斯泰勒法（Hilber-Hughes-Taylor method，HHT）［16-17］等。上述方法中僅最后2種為隱式積分法，其他均為顯式積分法。其中，對于可取較大步長值的仿真計算建議優(yōu)先選擇HHT和3 Gear法等隱式方法，而對于只能取小時間步長的仿真可考慮使用RK4和RKF45法等顯式方法。

此外，對于實(shí)際工程中的復(fù)雜剛?cè)狁詈舷到y(tǒng)，上述剛?cè)狁詈舷到y(tǒng)動力學(xué)方程的構(gòu)建方法往往難以應(yīng)用，因此，對于實(shí)際工程中復(fù)雜的剛?cè)狁詈蠁栴}，通常利用多剛體動力學(xué)和有限元法建立系統(tǒng)的剛?cè)狁詈夏Ｐ停M(jìn)而獲得系統(tǒng)的動力學(xué)響應(yīng)。

2? 基于區(qū)間K-L展開的剛?cè)狁詈舷到y(tǒng)動態(tài)不確定性序列模擬方法

2.1? 區(qū)間過程模型及區(qū)間K-L展開簡介

在傳統(tǒng)的隨機(jī)過程理論中，不確定過程X（t）在任意時刻t被視為一隨機(jī)變量，任意時刻變量之間的相關(guān)性由自相關(guān)函數(shù)來描述。然而隨機(jī)過程模型的構(gòu)建需要大量的樣本，在工程實(shí)際中很多時候難以獲取。而相對于精確的概率分布函數(shù)等統(tǒng)計信息，基于有限的樣本和工程經(jīng)驗(yàn)，參數(shù)的邊界信息往往相對容易獲得。如圖2所示，在區(qū)間過程模型中，任意時刻的不確定參數(shù)被描述成一個有界閉區(qū)間，不同時間點(diǎn)變量之間的相關(guān)性通過相應(yīng)的相關(guān)系數(shù)函數(shù)來描述［18-20］。

通常將區(qū)間過程{X（t）∈XI（t），t∈T}簡記為X（t）∈XI（t）或XI（t），其中T為參數(shù)t的參數(shù)集。關(guān)于參數(shù)t的函數(shù)XL（t）、XU（t）分別為區(qū)間過程X（t）∈XI（t）的下邊界（函數(shù)）和上邊界（函數(shù)），將區(qū)間過程XI（t）的樣本函數(shù)記為x（t）。區(qū)間過程X（t）∈XI（t）的中值函數(shù)和半徑函數(shù)分別定義為

Xm（t）=XU（t）+XL（t）2（11）

Xr（t）=XU（t）－XL（t）2（12）

對于區(qū)間過程XI（t），其中值函數(shù)Xm（t）反映該不確定過程隨時間t的整體變化趨勢。半徑函數(shù)Xr（t）非負(fù)，反映了XI（t）的波動程度。半徑函數(shù)越小，表明該區(qū)間過程的不確定性越弱，其中當(dāng)半徑函數(shù)Xr（t）=0時，則區(qū)間過程XI（t）退化為一確定性函數(shù)。XI（t）的協(xié)方差函數(shù)CovXIXI（ti，tj）和相關(guān)系數(shù)函數(shù)ρXIXI（ti，tj）=CovXIXI（ti，tj）Xr（ti）Xr（tj），ti，tj∈T用于表征區(qū)間過程在時間上的相關(guān)性。關(guān)于區(qū)間過程的詳細(xì)介紹可參考筆者前期研究工作［18-19］。

上述區(qū)間過程模型一般不能直接應(yīng)用于結(jié)構(gòu)不確定分析中，因而應(yīng)先將其通過數(shù)量有限的一組區(qū)間變量表示（即區(qū)間過程的離散）?；跁r間離散的方法為保證精度往往需要大量區(qū)間變量來表示時間的連續(xù)，但后續(xù)計算成本過高。除基于時間離散的方法以外，對于區(qū)間過程還可以通過級數(shù)展開的形式，將區(qū)間過程變換為無窮多個不相關(guān)的標(biāo)準(zhǔn)區(qū)間變量的確定性時間函數(shù)的組合，再保留有限重要項(xiàng)實(shí)現(xiàn)對區(qū)間過程的近似。

區(qū)間過程的區(qū)間K-L展開可表征為

X（t）=Xm（t）+∑∞j=1Xr（t）λjφj（t）ζj（13）

式中，ζj∈ζI=［－1，1］（j=1，2，…）為標(biāo)準(zhǔn)區(qū)間變量，且滿足∑∞i=1ζ2j≤1;λj∈［0，∞）、φj（t）分別為自相關(guān)系數(shù)函數(shù)ρXIXI（t，t′）的特征值與特征函數(shù)。

根據(jù)Mercer定理，對自相關(guān)系數(shù)函數(shù)ρXIXI（t，t′）有如下譜分解：

ρXIXI（t，t′）=∑∞j=1λjφj（t）φj（t′）（14）

其中，t，t′為任意的兩個時刻，特征值λj與特征函數(shù)φj（·）為如下第二類齊次方程Fredholm積分方程的解，即

∫TρXIXI（t，t′）φj（t）dt=λjφj（t′）（15）

正交特征函數(shù)滿足∫Tφi（t）φj（t）dt=δij，其中φi（t）為區(qū)別于φj（t）的特征函數(shù)，δij為Kronecker-delta函數(shù)。在式（15）中，令t=t′，此時有ρXIXI（t，t′）=1，等號兩邊在T上同時積分可得

∑∞j=1λj=|T|（16）

式中，|T|為時間參數(shù)集T的長度。

此外，區(qū)間過程的抽樣通過區(qū)間K-L展開實(shí)現(xiàn)。在實(shí)際工程中，利用無窮多個區(qū)間變量對隨時間變化的不確定參數(shù)進(jìn)行度量既不可能也無必要。在實(shí)際應(yīng)用時，可將式（13）所表示的級數(shù)項(xiàng)根據(jù)特征值的大小進(jìn)行降序排列，并對前M項(xiàng)進(jìn)行截斷。

記ζ（s）j（j=1，2，…，M）為標(biāo)準(zhǔn)區(qū)間變量ζj∈ζIj（j=1，2，…，M）的第s組樣本，滿足∑Mj=1ζ2j≤1，則區(qū)間過程X（t）∈XI（t）的樣本函數(shù)x（s）（t）可近似表示為

x（s）（t）=Xm（t）+∑Mj=1Xr（t）λjφj（t）ζ（s）j（17）

其中，特征值λj與特征函數(shù)φj（·）可根據(jù)式（15）確定。由于截斷近似帶來的誤差無法避免，故可將區(qū)間過程截斷近似的近似度κ定義為

κ=1|T|∑Mj=1λj（18）

其中，截斷項(xiàng)數(shù)M可根據(jù)近似度的要求進(jìn)行確定；為保證抽樣函數(shù)的準(zhǔn)確性，通常要求κ≥0.95。區(qū)間過程X（t）∈XI（t）的樣本函數(shù)x（s）（t）的具體抽樣步驟可參考文獻(xiàn)［18］。

2.2? 剛?cè)狁詈舷到y(tǒng)動態(tài)響應(yīng)上下邊界的序列模擬求解

在區(qū)間過程動態(tài)激勵作用下，剛?cè)狁詈舷到y(tǒng)的動力學(xué)響應(yīng)包括結(jié)構(gòu)振動和機(jī)構(gòu)運(yùn)動特性參數(shù)等，同樣具有上下邊界的動態(tài)不確定性［21］，如圖3所示。本文所研究的問題主要是對系統(tǒng)動態(tài)響應(yīng)上下邊界的求解。根據(jù)1.2節(jié)中所建立的剛?cè)狁詈蟿恿W(xué)方程（式（8）），下面結(jié)合區(qū)間K-L分解對輸入不確定激勵進(jìn)行表征，并給出該方程的非確定性描述形式。

首先，考慮由式（8）所確定的剛?cè)狁詈舷到y(tǒng)假設(shè)受X方向上區(qū)間過程激勵FX（t）。輸入激勵FX（t）由區(qū)間過程描述，即FX（t）∈FIX（t）=［FLX（t），F(xiàn)UX（t）］，F(xiàn)LX（t）、FUX（t）分別為不確定激勵FX（t）的下邊界函數(shù)和上邊界函數(shù)，F(xiàn)X（t）的自相關(guān)系數(shù)函數(shù)為ρXIXI（t，t′）。利用區(qū)間K-L展開對激勵FX（t）進(jìn)行表示，則有

FX（t）=FmX（t）+∑∞j=1FrX（t）λjφj（t）ζj（19）

將由式（19）所描述的激勵代入至剛?cè)狁詈蟿恿W(xué)方程（式（8））即可得到完整的非確定性描述形式。同理可得系統(tǒng)受到不同方向不確定激勵時的動力學(xué)方程，本文不再展開描述。

圖4展示出了用于尋找上邊界的序列模擬（sequential simulation，SS）方法的流程圖。序列模擬是通過一個序列過程來實(shí)現(xiàn)的，下邊界尋找的過程可通過類似的方式實(shí)現(xiàn)。在每個序列中，執(zhí)行直接蒙特卡羅（Monte Carlo，MC）模擬過程（即下述步驟（1）～步驟（2）），而每個序列中的采樣策略需要根據(jù)從先前序列獲得的結(jié)果響應(yīng)來更新［22］。序列模擬方法的詳細(xì)步驟如下：

（1）輸入?yún)^(qū)間過程X（t）∈XI（t）及其上下邊界函數(shù)、自相關(guān)系數(shù)函數(shù)，根據(jù)近似度的要求從M維單位球Ωζ={ζ|∑Mj=1ζ2j≤1}內(nèi)抽樣獲得標(biāo)準(zhǔn)區(qū)間變量ζj∈ζIj=［－1，1］的樣本ζ（s）j（j=1，2，…，M）；進(jìn)而通過區(qū)間K-L展開將標(biāo)準(zhǔn)區(qū)間變量ζj∈ζIj（j=1，2，…，M）的樣本ζ（s）j（j=1，2，…，M）映射為區(qū)間過程X（t）∈XI（t）的初始樣本函數(shù)x（s）（t）（s=1，2，…，ninitial），其中ninitial為初始樣本的個數(shù)。

（2）將生成的樣本輸入到所構(gòu)建的剛?cè)狁詈舷到y(tǒng)動力學(xué)方程或通過多剛體動力學(xué)和有限元法所建立的剛?cè)狁詈蟿恿W(xué)仿真模型中，求解得到輸出響應(yīng)為{y（x（s）（t），t），s=1，2，…，n}，n=nglobal+nlocal。nglobal、nlocal分別為上一個序列中新生成的全局樣本和局部樣本的數(shù)量。對于初始序列，有nglobal=ninitial，nlocal=0。

（3）從時間參數(shù)t的每個時刻{tj，tj∈［0，ts］}中求出所有樣本響應(yīng){y（x（s）（t），tj），s=1，2，…，n}對應(yīng)的最大值{y^U（tj），tj∈［0，ts］}，并獲得輸出響應(yīng)的上邊界y^U（t），其中ts為樣本響應(yīng)的總時間。

（4）識別對上邊界y^U（t）的貢獻(xiàn)因子，即識別在連續(xù)參數(shù)時間t上對整個時間歷程［0，ts］或部分時間間隔Δt具有貢獻(xiàn)最大值的響應(yīng)實(shí)現(xiàn){y（x（）（t），t）}，并對輸入的不確定參數(shù)尋找相應(yīng)的貢獻(xiàn)樣本{x（*）（t）}。

（5）如果從最近兩個序列獲得的上界之間的偏差足夠小，則認(rèn)為計算達(dá)到收斂并轉(zhuǎn)到步驟（6）；否則生成新的樣本（全局樣本和局部樣本），并轉(zhuǎn)到步驟（2）；全局樣本從整個有界空間XI（t）中產(chǎn)生，局部樣本在貢獻(xiàn)樣本的鄰域內(nèi)產(chǎn)生。

（6）將最近一次模擬序列中得到的上邊界作為近似解輸出。

上述方法為區(qū)間過程激勵下剛?cè)狁詈舷到y(tǒng)機(jī)構(gòu)運(yùn)動與結(jié)構(gòu)振動等動力學(xué)響應(yīng)邊界分析提供了一種通用的序列模擬策略。本節(jié)的后續(xù)內(nèi)容中，將討論序列模擬過程中的一些重要問題，如貢獻(xiàn)樣本的貢獻(xiàn)評估、貢獻(xiàn)樣本鄰域的定義、每個序列中局部樣本的數(shù)量、收斂標(biāo)準(zhǔn)的設(shè)置等。

在每個序列模擬中都需要識別出有貢獻(xiàn)的樣本。如前所述，通常大多數(shù)輸入樣本對輸出響應(yīng)的上邊界或下邊界沒有貢獻(xiàn)。例如，在每個序列模擬中獲得的上邊界必須由已有響應(yīng)中某些部分組成，因此，這些能夠構(gòu)成響應(yīng)邊界分量的輸入樣本就稱之為貢獻(xiàn)樣本。如圖5所示，將每一個貢獻(xiàn)樣本x（*）（t）的響應(yīng)實(shí)現(xiàn)記為y（x（）（t），t），其貢獻(xiàn)可以通過響應(yīng)實(shí)現(xiàn)覆蓋上邊界（或下邊界）的持續(xù)時間來評估，由此每個貢獻(xiàn)樣本對應(yīng)的貢獻(xiàn)指數(shù)定義為

γ［x（）（t）］=Δt［x（）（t）］ts（20）

其中，Δt［x（）（t）］為貢獻(xiàn)樣本x（）（t）的響應(yīng)實(shí)現(xiàn)覆蓋上邊界（或下邊界）的持續(xù)時間，圖5中x（1）（t）的響應(yīng)所覆蓋上邊界的時間非連續(xù)，即由兩個時間區(qū)域組成，“+”號表示兩個區(qū)域共同組成了該貢獻(xiàn)樣本響應(yīng)所覆蓋上邊界的總時間Δt［x（1）（t）］。

由序列模擬方法的過程可知，每個序列模擬中都將生成全局樣本和局部樣本（初始序列中沒有局部樣本）。全局樣本是在輸入?yún)^(qū)間不確定參數(shù)的整個有界空間內(nèi)生成的，而局部樣本則來自貢獻(xiàn)樣本的鄰域。全局樣本的數(shù)目可以根據(jù)要解決的數(shù)值問題的規(guī)模、非線性程度或計算花費(fèi)來給出。此外，在有貢獻(xiàn)的樣本附近的局部樣本被認(rèn)為更有可能產(chǎn)生響應(yīng)邊界值，因此，在下一個序列模擬中，優(yōu)選在貢獻(xiàn)樣本的鄰域內(nèi)采樣。

由上文可知，對于區(qū)間過程X（t）的抽樣即可等價于對有限個標(biāo)準(zhǔn)區(qū)間變量ζj∈ζIj=［－1，1］的抽樣；獲得區(qū)間變量的樣本ζ（s）j，借助于式（20）中的轉(zhuǎn)換即可近似獲得X（t）的樣本函數(shù)。同理可知，在貢獻(xiàn)樣本x（）（t）的鄰域內(nèi)進(jìn)行抽樣亦等價于在對應(yīng)的區(qū)間變量樣本ζ（s）j的鄰域內(nèi)進(jìn)行抽樣，進(jìn)而借助于式（17）中的轉(zhuǎn)換獲得局部樣本。圖6以二維標(biāo)準(zhǔn)區(qū)間變量空間為例對不同形狀的貢獻(xiàn)樣本鄰域進(jìn)行介紹。圖6a所示的圓形鄰域是一個圓形區(qū)域，對于多個標(biāo)準(zhǔn)隨機(jī)變量，鄰域變?yōu)橐载暙I(xiàn)樣本為中心、以r*為半徑的多維球形區(qū)域。圖6b所示的正方形鄰域，它為以貢獻(xiàn)樣本為中心、邊長為2r*的正方形或多維立方體形區(qū)域。為了便于描述，本文引入比率β來表示鄰域相對于整個有界空間的大小。對于圓形鄰域和方形鄰域，定義為 r*=βMV（Ωζ），其中V（Ωζ）為有界空間的體積或面積；w（ζI）為區(qū)間變量的長度，對于標(biāo)準(zhǔn)區(qū)間變量，有w（ζI）=2。通過測試一些區(qū)間不確定性問題，

本文建議將比率β設(shè)置在1/12到1/4之間。

各個貢獻(xiàn)樣本鄰域內(nèi)的局部樣本數(shù)nlocal［x（）（t）］可根據(jù)貢獻(xiàn)指數(shù)γ［x（）（t）］給出，且有nlocal［x（）（t）］∝γ［x（）（t）］。

定義y-（k）（t）為第k次序列模擬所獲得的上邊界，如圖7所示，第k次模擬序列的絕對誤差可表示為ε（k）（t）=|y-（k）（t）－y-（k－1）（t）|。則收斂準(zhǔn)則可以描述為：

maxt ε（k）（t）≤ε0（21）

式中，ε0為用于控制誤差的預(yù)設(shè)值。

3? 算例分析

3.1? 振動基柔性臂

考慮圖8所示的振動基柔性臂系統(tǒng)，坐標(biāo)系OXY表示系統(tǒng)的慣性坐標(biāo)系，在柔性臂上建立浮動坐標(biāo)系oxy，其中柔性臂可以o點(diǎn)為旋轉(zhuǎn)中心翻轉(zhuǎn)至指定角度，柔性臂的末端與電機(jī)軸固接，屬于懸臂梁模型約束條件，θ為翻轉(zhuǎn)部分相對于X軸的夾角。x為未變形時柔性臂的軸線。X軸所在平面表示系統(tǒng)未發(fā)生垂直振動時所在水平面，O點(diǎn)處為基座質(zhì)心；Y軸為基座質(zhì)心的垂直振動的方向?；馁|(zhì)量取mz=0.5 kg，柔性臂的參數(shù)取長度L=8 m，截面積S=7.3×10－3? m2，慣性矩I=8.22×10－9 m4，

密度ρ=2.7667×103 kg/m3，彈性模量E=68.952 GPa。

本算例同時考慮柔性臂結(jié)構(gòu)柔性和基座垂直振動，其中，基座底部受到一動態(tài)載荷f（t），且被度量為一中值為0、半徑為2 N的平穩(wěn)區(qū)間過程f（t）∈fI（t），其相關(guān)性系數(shù)函數(shù)為ρXIXI（t，t′）=exp（－|t－t′|/lt），其中l(wèi)t為相關(guān)長度，本文取lt=0.3 s。對于指數(shù)型自相關(guān)系數(shù)函數(shù)，相關(guān)長度lt越大，區(qū)間過程f（t）的時間相關(guān)性越強(qiáng)，樣本函數(shù)相對也越平滑［15］。與此同時柔性臂受一確定性力矩驅(qū)動下繞o點(diǎn)旋轉(zhuǎn)，且驅(qū)動力矩為

τ（t）=τ0sin2πtTm? 0≤t≤Tm

0t>Tm（22）

式中，Tm為力矩作用時間，取Tm=10 s；τ0為力矩幅值，取τ0=5 N·m。

在實(shí)際中，模態(tài)坐標(biāo)一般由前幾階模態(tài)主導(dǎo)，本算例為單自由度柔性臂的振動模型，取前兩階模態(tài)即可得到精準(zhǔn)的結(jié)果［9］。且假定系統(tǒng)只在Y方向做平動，采取廣義坐標(biāo)q=（Y，θ，A，B），其中，A（t）=（A1（t），A2（t））T，B（t）=（B1（t），B2（t））T。

使用第1節(jié)所述剛?cè)狁詈舷到y(tǒng)動力學(xué)建模方法，可得到如下方程：

M22M23M24M25

M32M33M34M35

M42M43M440

M52M530M55

Y¨θ¨A¨B¨=QYQθQAQB（23）

本算例中的柔性臂模態(tài)取前兩階，則其橫向變形ω2（x，t）可描述為如下模態(tài)方程：

ω2（x，t）=∑2i=1Φ（i）y（x）Bi（t）=

Φ（1）y（x）B1（t）+Φ（2）y（x）B2（t）（24）

基于上述剛?cè)狁詈戏匠探M，本算例分別通過直接蒙特卡羅模擬方法（以下簡稱MC法）和序列模擬方法（以下簡稱SS法）分別計算了柔性臂末端橫向位移ω2與柔性臂翻轉(zhuǎn)角θ隨時間變化的上下邊界，將仿真時間設(shè)為t=15 s，且仿真的初始條件均為0。為了保證結(jié)果的精度，本文將MC法的模擬次數(shù)設(shè)置為20 000。在調(diào)用SS法對該問題進(jìn)行求解時，響應(yīng)上下邊界的收斂的最大誤差均設(shè)置為maxt ε（k）（t）≤ε0=1×10－7，且應(yīng)用比率β=0.2的多維球形作為鄰域。應(yīng)用MC法與SS法對上述問題進(jìn)行求解，所得結(jié)果如圖9和圖10所示，結(jié)果表明MC法所得邊界與SS法所得上下邊界高度重合，甚至比MC法所得邊界更寬。與基于區(qū)間運(yùn)算的區(qū)間不確定性分析方法不同，SS法和MC法均為非嵌入式方法，且給出近似內(nèi)部解，得到的響應(yīng)邊界將嚴(yán)格地被包絡(luò)于實(shí)際結(jié)果之內(nèi)，因此，計算的響應(yīng)邊界越寬，表明越接近精確響應(yīng)邊界。表1、表2所示為不同時刻下SS法所求邊界相比于MC法所求邊界的誤差，結(jié)果表明不同時刻下SS法相對于MC法的誤差最大為7.69%。與此同時，相比于MC法中對原模型進(jìn)行了20 000次調(diào)用，SS法在求解上下邊界過程中分別對原剛?cè)狁詈戏匠踢M(jìn)行了380次和222次的調(diào)用，分別僅有MC法調(diào)用次數(shù)的1/50和1/100左右。由上述討論可知，在該剛?cè)狁詈蠁?/p>

題中，SS法能夠以較小的計算量獲得較高精度的解。

3.2? 剛?cè)狁詈线\(yùn)動機(jī)構(gòu)

依托某工程實(shí)際剛?cè)狁詈线\(yùn)動機(jī)構(gòu)，建立圖11所示剛?cè)狁詈线\(yùn)動機(jī)構(gòu)模型，其中，桿1～桿3和質(zhì)量塊A、B、C、D均為剛體，桿4為柔性體。坐標(biāo)系OXYZ表示系統(tǒng)的慣性坐標(biāo)系，在桿4上建立浮動坐標(biāo)系oxy。桿1和桿2、桿3和桿4分別在B、D處固連；桿2在C處與桿3鉸接且可繞質(zhì)量塊C的質(zhì)心C點(diǎn)旋轉(zhuǎn)，轉(zhuǎn)角θ為桿3轉(zhuǎn)動后與OY軸的夾角。桿1～桿4的長度依次為40 mm、60 mm、70 mm、120 mm，且直徑均為6 mm；質(zhì)量塊A的尺寸為10 mm×10 mm×20 mm；其余質(zhì)量塊尺寸均為10 mm×10 mm×10 mm。假定質(zhì)量塊A底部中心受到一不確定性激勵f（t），且在其作用下機(jī)構(gòu)產(chǎn)生的時變加速度a可通過底部加速度傳感器進(jìn)行測量，并以區(qū)間過程模型進(jìn)行度量，即a（t）∈aI（t）=［－2，2］m/s2。且桿3繞C點(diǎn)處轉(zhuǎn)動的運(yùn)動規(guī)律為θ=π6sin（20πt）rad。本算例仿真分析選用的材料參數(shù)列入表3。

使用UG和ABAQUS軟件建立剛?cè)狁詈夏Ｐ蜁r，桿4是重點(diǎn)關(guān)注的部分，作柔性化處理生成柔性體。其余部分作為剛性體建模，然后將剛性和柔性部分裝配好，形成剛?cè)狁詈夏Ｐ停鐖D12所示。

分別使用蒙特卡羅直接模擬程序和前文所述的序列模擬程序并聯(lián)合ADAMS仿真軟件對系統(tǒng)動力學(xué)進(jìn)行分析，仿真時間t=5 s。為了保證結(jié)果的精度，本文將MC法的模擬次數(shù)設(shè)置為10 000。調(diào)用SS方法對該問題進(jìn)行求解時，響應(yīng)上下邊界的收斂的最大誤差設(shè)置為ε0=1×10－3 mm，且采用比率β=0.15的多維立方體形鄰域；SS法的模擬總數(shù)為201次。比較了SS法和MC法評估桿4的末端橫向位移ω2響應(yīng)圖的上邊界和下邊界，仿真結(jié)果如圖13所示。其中，不同時刻下SS法所求邊界相比于MC法所求邊界的誤差如表4所示。

the end of a flexible rod

如圖13和表4所示，雖然SS法模擬次數(shù)要少得多，但是SS法仍給出了更好的解。在這種情況下，SS法給出的上限甚至比MC法給出的上限更高，下限比MC法給出的下限更低。如前所述，MC法和SS法均可獲得響應(yīng)近似邊界，所得響應(yīng)界限將嚴(yán)格包含在實(shí)際結(jié)果內(nèi)，且表明響應(yīng)界限越寬，越接近準(zhǔn)確結(jié)果。由此可以看出，SS法比MC法使用的模擬次數(shù)要少得多，且結(jié)果更好，故SS法在該數(shù)值示例中顯示出了更高的計算效率。

3.3? 某型碼垛機(jī)器人執(zhí)行機(jī)構(gòu)

ABB的某型號碼垛機(jī)器人的最大載重量為450 kg，作業(yè)半徑為3.2 m，該本體結(jié)構(gòu)見圖14。其結(jié)構(gòu)為典型的雙平行四邊形結(jié)構(gòu)，主要由基座、腰部、腕部和臂部機(jī)構(gòu)組成［23］，臂部機(jī)構(gòu)按照一定比例組成平行四邊形，具有4個相互獨(dú)立的自

由度，分別為腰部相對于底座的回轉(zhuǎn)運(yùn)動、大臂的前后仰俯運(yùn)動、小臂的上下仰俯運(yùn)動、腕部的回轉(zhuǎn)運(yùn)動。

碼垛機(jī)器人臂部結(jié)構(gòu)為懸臂的細(xì)長桿件，在極端振動工況下的運(yùn)動過程中會出現(xiàn)相對較大的柔性變形。由于碼垛機(jī)器人臂部柔性變形對機(jī)器人的運(yùn)行精度及穩(wěn)定性影響較大，以剛性件對其進(jìn)行分析難以滿足使用要求，因此采用剛?cè)狁詈戏椒筛訙?zhǔn)確地分析碼垛機(jī)器人的運(yùn)行狀態(tài)。

使用UG和ABAQUS軟件建立剛?cè)狁詈夏Ｐ蜁r，小臂是重點(diǎn)關(guān)注的部分，作柔性化處理生成柔性體。其余部分作為剛性體建模，然后將剛性和柔性部分裝配好，形成剛?cè)狁詈夏Ｐ?。如圖15所示，忽略螺釘、墊片、軸套等對系統(tǒng)影響較小的零件，簡化各軸、電機(jī)、減速機(jī)等，然后設(shè)置單位、重力加速度和材料屬性等參數(shù)。

考慮機(jī)構(gòu)工作時的極端振動工況，假定底座中心受到一垂直不確定性激勵f（t）。為模擬此振動工況，本算例以加速度a（t） 的形式對模型進(jìn)行驅(qū)動，且被度量為一中值為0、半徑為20 m/s2的平穩(wěn)區(qū)間過程，其相關(guān)性系數(shù)函數(shù)為ρXIXI（t，t′）=exp（－|t－t′|/lt），取lt=0.4 s。機(jī)構(gòu)一個工作周期的工作軌跡如圖16所示：從0～0.5 s，機(jī)器人末端執(zhí)行器向下移動700 mm，然后停頓0.5 s抓取400 kg的載荷；從1.0～1.5 s，末端執(zhí)行器向上抬高700 mm；然后腰部旋轉(zhuǎn)90°同時伴隨手臂運(yùn)動達(dá)到碼放位置，這些動作在1.5～3.0 s之間完成；之后用0.5 s碼放負(fù)載；最后用0.5 s返回到初始位置。完成一個完整的作業(yè)總共用時4 s。

分別使用MC法和上文所述的SS法并聯(lián)合ADAMS仿真軟件對上述仿真條件下機(jī)器人末端執(zhí)行器的速度、加速度曲線進(jìn)行分析，仿真時間t=4 s。

基于上述剛?cè)狁詈舷到y(tǒng)仿真分析，本算例分別通過MC法和SS法依次計算了機(jī)器人末端執(zhí)行器的速度、加速度隨時間變化的上下邊界。為了保證結(jié)果的精度，本文將MC法的模擬次數(shù)設(shè)置為10 000。調(diào)用SS法對該問題進(jìn)行求解時，響應(yīng)上下邊界的收斂的最大誤差均設(shè)置為ε0=1×10－4，且應(yīng)用比率β=0.15的多維球形作為鄰域。應(yīng)用MC法和SS法對上述問題進(jìn)行求解，所得結(jié)果如圖17和圖18所示，結(jié)果表明SS法所得邊界與MC法所得上下邊界高度重合。表5和表6所示為不同時刻下SS法所求邊界相比于MC法所求邊界的誤差，結(jié)果表明兩者最大誤差為9.76%。與此同時，相比于MC法中對系統(tǒng)仿真模型進(jìn)行了10 000次調(diào)用，SS法在求解上下邊界過程中分別對原剛?cè)狁詈戏匠踢M(jìn)行了333次和376次的調(diào)用，分別僅有MC法調(diào)用次數(shù)的3%與5%左右。

如圖17、圖18和表5、表6所示：在整個工作周期內(nèi)，加速度波動較大，說明機(jī)器人在運(yùn)動過程中存在較大振動，尤其是在運(yùn)動狀態(tài)發(fā)生改變時，運(yùn)動極不平穩(wěn)，從而導(dǎo)致機(jī)器人運(yùn)動精度降低。且在不同時刻下SS法所得邊界與MC法模擬

10 000次后的邊界高度重合甚至更寬。如前所述，這兩種方法計算所得均為近似內(nèi)部解，得到的響應(yīng)邊界將嚴(yán)格地被包絡(luò)于實(shí)際結(jié)果之內(nèi)，計算所得響應(yīng)邊界越寬，表明越接近精確響應(yīng)邊界，因此在較為復(fù)雜的剛?cè)狁詈蠁栴}中，序列模擬方法仍能夠以較小的計算量獲得較高精度的解。

4? 結(jié)論

本文針對動態(tài)不確定性激勵作用下的剛?cè)狁詈舷到y(tǒng)動力學(xué)問題提出了一種基于區(qū)間過程模型的動力學(xué)不確定性序列模擬方法，可高效計算結(jié)構(gòu)振動與機(jī)構(gòu)運(yùn)動耦合動力學(xué)響應(yīng)的上下邊界?？紤]到直接蒙特卡羅模擬方法在進(jìn)行系統(tǒng)動力學(xué)響應(yīng)上下邊界時，因絕大多數(shù)區(qū)間過程抽樣及基于此的剛?cè)狁詈蟿恿W(xué)分析模擬對機(jī)構(gòu)運(yùn)動狀態(tài)參數(shù)和結(jié)構(gòu)振動響應(yīng)的邊界求解不具有任何作用，導(dǎo)致計算效率極低的問題，本文方法利用序列抽樣模擬的方式逐步逼近系統(tǒng)在整個時域上的動態(tài)響應(yīng)上下邊界。該方法通過識別當(dāng)前剛?cè)狁詈?/p>

動力學(xué)計算模擬序列中對動態(tài)響應(yīng)上邊界或下邊界具有貢獻(xiàn)的區(qū)間過程參數(shù)樣本集，作為下一模擬序列中的局部加密抽樣中心，從而有效提高了響應(yīng)邊界的計算收斂性。算例驗(yàn)證結(jié)果表明，對于剛?cè)狁詈舷到y(tǒng)大范圍運(yùn)動及振動響應(yīng)上下邊界的求解，相比于直接蒙特卡羅模擬方法，序列模擬方法能以較小的計算量獲得較高精度的解。

參考文獻(xiàn)：

［1］? 陳思佳，剛?cè)狁詈蠁栴}與空間多桿柔性機(jī)械臂的動力學(xué)建模理論研究［D］. 南京：南京理工大學(xué)，2012.

CHEN Sijia.Research on Dynamic Modeling Theory of Rigid-flexible Coupling Problem and Spatial Multi-rod Flexible Manipulator［D］. Nanjing：Nanjing University of Science and Technology， 2012.

［2］? 王斌銳， 方水光， 嚴(yán)冬明. 機(jī)器人手臂的剛?cè)狁詈辖＜皵[動模態(tài)對比［J］.中國機(jī)械工程，2012，23（17）：2092-2097.

WANG Binrui， FANG Shuiguang， YAN Dongming. Rigid-flex Coupling Modeling and Swing Mode Comparison of Robot Arm［J］.China Mechanical Engineering，2012，23（17）：2092-2097.

［3］? 莊未， 方琛瑋， 劉曉平. 1P5R柔性關(guān)節(jié)機(jī)械臂動力學(xué)分析與控制［J］.中國機(jī)械工程，2009，20（24）：2907-2911.

ZHUANG Wei， FANG Chenwei， LIU Xiaoping. Dynamic Analysis and Control of 1P5R Flexible Joint Manipulator［J］.China Mechanical Engineering， 2009， 20（24）：2907-2911.

［4］? 吳勝寶，章定國. 大范圍運(yùn)動剛體-柔性梁剛?cè)狁詈蟿恿W(xué)分析［J］. 振動工程學(xué)報，2011，24（1）：1-7.

WU Shengbao， ZHANG Dingguo. Dynamic Analysis of Rigid-flexible Coupling of Large-scale Motion Rigid Body-flexible Beam［J］. Journal of Vibration Engineering， 2011， 24（1）：1-7.

［5］? 胡啟國，周松，考慮剛?cè)狁詈系墓I(yè)機(jī)器人多目標(biāo)可靠性拓?fù)鋬?yōu)化［J］. 計算機(jī)集成制造系統(tǒng)，2020，26（3）：623-631.

HU Qiguo， ZHOU Song. Multi-objective Reliability Topology Optimization of Industrial Robots Considering Rigid-flexible Coupling［J］. Computer Integrated Manufacturing Systems， 2020，26（3）：623-631.

［6］? 辛鵬飛，榮吉利，項(xiàng)陽，等. 柔性空間機(jī)械臂區(qū)間參數(shù)不確定性分析［J］. 北京理工大學(xué)學(xué)報，2017，37（10）：1056-1060.

XIN Pengfei， RONG Jili， XIANG Yang， et al. Uncertainty Analysis of Interval Parameters of Flexible Space Manipulator［J］. Journal of Beijing Institute of Technology， 2017，37（10）：1056-1060.

［7］? WANG L， YANG G， XIAO H， et al. Interval Optimization for Structural Dynamic Responses of an Artillery System under Uncertainty［J］. Engineering Optimization， 2020， 52（2）：343-366.

［8］? 李勇. 某供輸彈機(jī)供藥機(jī)構(gòu)的可靠性研究［D］.南京：南京理工大學(xué)，2017.

LI Yong. Research on the Reliability of the Medicine Supply Mechanism of a Certain Bomb Feeding Machine［D］. Nanjing：Nanjing University of Science and Technology， 2017.

［9］? 宋清林.考慮間隙和柔性體的火炮協(xié)調(diào)器動力學(xué)與可靠性分析［D］.沈陽：東北大學(xué)，2018.

SONG Qinglin. Dynamics and Reliability Analysis of Artillery Coordinator Considering Gaps and Flexible Bodies［D］. Shenyang：Northeastern University， 2018.

［10］? YANG C， LU W， XIA Y. Reliability-constrained Optimal Attitude-vibration Control for Rigid-flexible Coupling Satellite Using Interval Dimension-wise Analysis［J］. Reliability Engineering & System Safety， 2023， 237：109382.

［11］? 郝秀清，胡福生，陳建濤.基于牛頓歐拉法的3PTT并聯(lián)機(jī)構(gòu)動力學(xué)分析及仿真［J］.中國機(jī)械工程，2006（增刊2）：32-36.

HAO Xiuqing， HU Fusheng， CHEN Jiantao. Dynamic Analysis and Simulation of 3PTT Parallel Mechanism Based on Newton-Euler Method［J］.China Mechanical Engineering，2006（S2）：32-36.

［12］? 朱暢. 基于奇異攝動法的振動基柔性機(jī)械臂的非線性魯棒控制研究［D］. 武漢：武漢科技大學(xué)，2018.

ZHU Chang. Research on Nonlinear Robust Control of Vibration-based Flexible Manipulator Based on Singular Perturbation Method［D］. Wuhan：Wuhan University of Science and Technology， 2018.

［13］? 郭晛. 受約束柔性多體系統(tǒng)動力學(xué)數(shù)值計算方法若干問題研究［D］. 南京：南京理工大學(xué)，2019.

GUO Xian. Research on Several Issues in Numerical Calculation Methods for Constrained Flexible Multi-body System Dynamics［D］. Nanjing：Nanjing University of Science and Technology， 2019.

［14］? 袁兆鼎，費(fèi)景高，劉德貴. 剛性微分方程初值問題數(shù)值解法［M］. 北京：科學(xué)出版社，1987.

YUAN Zhaoding， FEI Jinggao， LIU Degui. Numerical Solution to Initial Value Problem of Rigid Differential Equations［M］. Beijing：Science Press， 1987.

［15］? NEWMARK N M. A Method of Computation for Structural Dynamics［J］. Journal of the Engineering Mechanics Division，1959，85（3）：67-94.

［16］? HILBER H， HUGHES T， TAYLOR R. Improved Numerical Dissipation for Time Integration Algorithms in Structural Dynamics［ J ］. Earthquake Engineering & Structural Dynamics，1977，5（2）：283-292.

［17］? NEGRUT D， RAMPALLI R， OTTARSSON G， et al. On an Implementation of the Hilber-Hughes-Taylor Method in the Context of Index 3 Differential-algebraic Equations of Multibody Dynamics［J］. Journal of Computational and Nonlinear Dynamics， 2007，2（1）：73-85.

［18］? 倪冰雨. 區(qū)間過程與區(qū)間場模型及在結(jié)構(gòu)不確定性分析中的應(yīng)用［D］.長沙：湖南大學(xué)，2018.

NI Bingyu. Interval Process and Interval Field Models and Their Application in Structural Uncertainty Analysis［D］. Changsha：Hunan University， 2018.

［19］? JIANG C， NI B， HAN X， et al．Non-probabilistic Convex Model Process：a New Method of Time-variant Uncertainty Analysis and Its Application to Structural Dynamic Reliability Problems， Computer Methods in Applied Mechanics and Engineering， 268 （2014） 656-676.

［20］? JIANG C， HAN X， LU G Y， et al． Correlation Analysis of Non-probabilistic Convex Model and Corresponding Structural Reliability Technique， Computer Methods in Applied Mechanics and Engineering， 2011，200：2528-2546.

［21］? 曾光， 段民封， 倪冰雨， 等. 機(jī)械臂架系統(tǒng)的非隨機(jī)振動分析［J］. 中國機(jī)械工程， 2019， 30（10）：1142.

ZENG Guang， DUAN Minfeng， NI Bingyu， et al. Non-random Vibration Analysis of Robotic Boom System［J］. China Mechanical Engineering， 2019， 30（10）：1142.

［22］? NI B Y， JIANG C， WU P G， et al. A Sequential Simulation Strategy for Response Bounds Analysis of Structures with Interval Uncertainties［J］. Computers & Structures， 2022， 266：106785.

［23］? 陳亞梅. 基于ADAMS的碼垛機(jī)器人動力學(xué)剛?cè)狁詈戏治觯跩］. 包裝工程，2021，42（17）：261-265.

CHEN Yamei. Rigid-flexible Coupling Analysis of Palletizing Robot Dynamics Based on ADAMS［J］. Packaging Engineering， 2021， 42（17）：261-265.

（編輯? 胡佳慧）

作者簡介：

劉延浩，男，1999年生，碩士研究生。研究方向?yàn)闄C(jī)械結(jié)構(gòu)動力學(xué)分析。E-mail：lyh1@hnu.edu.cn。

倪冰雨（通信作者），男，1989年生，副教授。研究方向?yàn)榻Y(jié)構(gòu)不確定性分析。E-mail：nby@hnu.edu.cn。