剛柔耦合接觸動態粘合算法研究

高 浩，干 鋒，戴煥云

(西南交通大學 牽引動力國家重點實驗室，成都 610031)

剛柔耦合多體系統中包含了剛體、彈性體以及體與體之間的連接力元和鉸，剛柔耦合仿真技術已經基本成熟，幾乎所有通用多體動力學軟件都支持剛柔耦合仿真，一般采用模態疊加法近似描述彈性體的動態響應。模態疊加法是基于小變形假設的線性處理方法，忽略了結構的非線性特性，同時模態截斷造成了高頻振動信息的丟失。對于非線性彈性變形或局部模態對仿真結果影響很大的系統，彈性體需要直接采用有限元瞬態動力學分析?？梢酝ㄟ^一定的算法將多系統中多剛體和彈性體分為不同的子系統，分別在不同的軟件中進行仿真。這是一種松耦合方式，關鍵技術是如何通過數學工具建立嚴密的接口關系，通過接口傳遞界面力和運動信息包括位移、速度和加速度，并保證接口界面處同時滿足平衡條件和協調條件。剛柔耦合接觸問題是剛柔耦合系統動力學仿真的難題之一，同時也有很多工程應用如齒輪接觸動力學和車輛系統弓網接觸動力學研究。文獻[4]中采用有限元預計算的方法考慮齒輪的彈性變形，研究剛柔耦合特性對齒輪嚙合傳動特性的影響規律。文獻[5]采用直接有限元積分的方法對弓網接觸動力學進行研究。Wang[1－2]提出了一種剛柔耦合仿真的松耦合方法稱為粘合算法，子系統界面之間傳遞的信息通過協調器即粘合算法處理并傳遞給與之相連的其他子系統，這種方法將子系統的運動信息作為輸入，輸出的界面力作為子系統的下一步積分的邊界條件。該方法可以用于多體系統剛柔耦合仿真分析。Run[3]在此基礎上改進了粘合算法，針對多體動力學中彈性鉸節點提出了彈性粘合算法，減小了連接界面的數值振蕩，提高了算法的穩定性。但是粘合算法是針對連接界面固定的子系統提出的，不適合連接界面變化接觸問題。為此本文在粘合算法的基礎上提出了動態粘合算法，即粘合算法中的裝配矩陣和柔度矩陣是根據接觸位置動態變化的，以保證接觸界面平衡條件和協調條件同時滿足。并將其用于處理剛體與彈性體或彈性體與彈性體之間的動態接觸問題。當接觸發生在沒有節點的位置時，如何處理不匹配的網格非常關鍵，一般采用形函數插值法[6]，Kim[7－8]提出界面單元法，將無網格法中的移動最小二乘近似引入到有限元中。本文在粘合矩陣中考慮單元形函數使得動態粘合算法可以處理非節點位置接觸的情況。

1 多體系統粘合算法[2]

多體系統一般有剛體或彈性體多個部件通過力元、鉸接或約束組合而成，假設部件i和部件i+1通過界面節點p和p'連接，如圖1所示。

圖1 多體系統部件粘合Fig.1 Components gluing of multibody system

要保證平衡條件和協調條件同時滿足，首先假設合適的界面力向量，使得每個子系統滿足平衡條件，在仿真過程中只需要考慮協調條件即界面連接點的運行信息，包括位移、速度和加速度。界面力是運動信息的函數，如式所示:

式(1)可表示為:

接觸界面應當滿足協調條件，界面運動信息裝配向量為U，則:

其中表示一對連接點的位移，假設e為協調條件的誤差向量，界面應滿足協調條件:

求解式(3)中的界面力向量F，使其滿足協調條件式(5)，利用牛頓法－拉夫遜法求解該非線性方程組式，迭代關系式為:

其中:k為迭代次數，Λ稱為粘合矩陣，Fk和Fk+1分別為第k步和第k+1步的界面力，具體計算方法如下:

其中:B和C稱為裝配矩陣，其元素由0、－1和1組成;G稱為界面柔度矩陣，式(10)中m和n分別為G的第m行，第n列元素，利用式(1)和式(10)計算獲得。對于線性系統可近似的通過在連接界面各自由度施加單位載荷獲取柔度矩陣，則其中為在連接界面自由度m施加單位載荷時產生的位移為不加任何載荷時的初始位移。i－th表示i部件對應的元素。

通過式(6)計算獲得的界面力，以外載荷的形式施加到各部件的動力學方程中，實現各部件之間的粘合。

2 接觸動態粘合法則

兩個部件i和i+1相互接觸，其中i為接觸體，i+1為目標體，假設i結構上接觸點為p點，對應的部件i+1上的p'位置，p與p'為接觸對，如圖2所示。若p'不是節點，p'坐標由該點所在單元的節點坐標和形函數插值得到，假設接觸單元為n個節點單元，p'點坐標如式(11)所示:

其中:Ni接觸單元中節點i的形函數，ui為節點i的坐標，Xp'為p'在接觸單元中的自然坐標向量。利用粘合算法求解兩個部件之間的接觸問題，首先要確定連接界面，接觸體上的接觸點p和目標體中對應點p'所在單元相關的節點組成連接界面，兩個部件之間的通過非線性接觸力元連接，接觸剛度為Kc，隨著接觸狀態的變化，接觸界面是動態變化的，因此本文發展動態粘合算法來求解接觸問題。

接觸點之間不是固接的，令接觸點之間的距離qc為:

其中界面裝配誤差向量為:

圖2 兩接觸部件粘合Fig2.Gluing of two contact components

協調方程可寫成:

假設接觸點之間為非線性接觸彈簧，接觸力為:

式(7)中的粘合矩陣改寫為:

動態粘合算法需要動態的更新裝配矩陣和界面柔度矩陣，裝配矩陣B和C的元素不再只有1、－1和0，而是和單元的形函數相關，如式(17)所示。

其中:N=[N1(Xp')，N2(Xp')，…，Nn(Xp')]，接觸界面柔度矩陣G也需要動態更新，將式(11)代入式(1)中，結合式(10)求得G矩陣。在整個接觸過程中裝配矩陣和界面柔度矩陣的維數和數值都會發生變化，因此每一積分步都需要對這些矩陣進行更新。

3 實例仿真

使用本文發展的動態粘合算法，對如圖3所示的系統進行動力學仿真，該系統包含簡支梁和彈簧－質量塊兩個部件，部件之間通過接觸單元連接，mc為質量塊的質量，yc為質量塊運動自由度，ks為支撐剛度，kc為接觸剛度，ys為支撐基礎的運動，V0為彈簧－質量塊向前運動速度。

選取仿真參數如表1所示。

圖3 簡支梁和質量－彈簧系統接觸模型Fig3.Simple-supported beam and mass-spring systems contact model

表1 動力學仿真參數Tab.1 Parameters of dynamic simulation

支撐基礎正弦輸入為 ys=0.001+0.002sin(4πt)，在該正弦輸入激擾下，簡支梁和質量塊會出現接觸和脫開的狀態轉換。采用統一方程法和動態粘合算法兩種方法求解該系統的響應。

3.1 統一方程法

采用有限元法建立平面梁的動力學方程，梁單元有兩個節點，每個節點有水平、垂向和轉動三個自由度，選用梁單元的形函數為:

其中:l為單元長度，ξ為梁單元坐標系下的坐標。整個系統的統一動力學方程為:

N下標表示式(18)形函數中對應位置的元素，ξp'表示p'點在梁單元坐標系中坐標。NI的非零元素分別對應節點 i和節點 j的垂向和轉動自由度，即 uiy，uiθ，ujy和ujθ。用Newmarker法求解式(19)，每個積分步判斷接觸狀態，當接觸脫開時令kc=0。

3.2 動態粘合算法

系統中包含兩個部件，平面簡支梁和彈簧－質量系統，兩個部件只有一個連接界面，即接觸對，兩個部件的動力學方程分別為:

假設彈簧質量系統為部件1，簡支梁為部件2，部件1的界面位移和界面力為，部件 2 的界面位移和界面力為[FiyMiθFjyMjθ]T，其中 uiy，uiθ，ujy和 ujθ分別為接觸單元i節點和 j節點的垂向和轉動位移，Fiy，Miθ，Fjy和Mjθ分別為作用在接觸單元i節點和j節點的垂向節點力和扭矩。

兩個部件之間存在一個連接界面，因此界面裝配協調條件為一個方程:

式中:N’(ξp')=[N22(ξp') N23(ξp') N25(ξp')N26(ξp')]。接觸點之間的距離為q=up－up'，此時協調方程為e－q=0，接觸力為FcA=kcq，根據有限元法將集中力轉化成等效節點力，即 FcR=－N’(ξp')FcA，取式(4)中的F=FcA，根據式(8)、式(9)和式 (17)計算部件1和部件2的裝配矩陣為:

式(24)～式(27)代入式(16)中，粘合矩陣可表示為:

界面力迭代關系式為:

4 仿真結果及討論

用統一方程法和動態粘合算法對上述系統進行動力學仿真，圖4為質量塊運動軌跡及簡支梁在不同時刻的變形。質量塊的垂向位移時程曲線如圖5，梁中部節點的垂向位移時程曲線如圖6所示，兩種方法計算位移結果基本一致。圖7為兩種算法計算的位移結果的絕對誤差，可見誤差非常小，最大誤差發生在由脫離狀態到接觸狀態變化瞬間，穩定接觸過程中絕對誤差小于5e－7 m，完全可以滿足仿真精度需求。圖8為梁和質量塊之間的接觸力，兩種算法在接觸力突變時仍保持很好的一致性，說明粘合算法比較穩定。

圖4 不同時刻接觸系統的狀態Fig.4 Status of contact system at different time

圖5 質量塊垂向位移Fig.5 Vertical displacement of mass block

圖6 梁中部垂向位移Fig.6 Vertical displacement of central node of the beam

圖7 位移絕對誤差Fig.7 Absolute error of displacement

圖8 梁與質量塊之間的接觸力Fig.8 Contact force between beam and mass block

在使用粘合算法時，用上一積分步的界面力作為界面力迭代的初始條件，經過1～2次迭代即可滿足仿真精度要求?？梢妱討B粘合算法可以很好的替代統一方程求解，可以將復雜系統分解成多部件，分別求解在不同的軟件或計算機上進行求解，實現分布式計算。

對比式(19)和(22)可以發現接觸問題為狀態非線性問題，使用統一方程求解時，剛度矩陣是變化的需要動態更新，無法進行事先分解。使用動態粘合算法，剛度矩陣不變，對以剛度矩陣規模較大的有限元系統，可進行事先分解，大大提高計算速度。

5 結論

(1)本文發展了一種動態粘合算法，可以用來解決剛性體與彈性體或彈性體與彈性體之間的接觸問題，該方法通過動態的構建粘合矩陣，并迭代求解一組界面力，使得兩個接觸的系統在接觸界面同時滿足平衡條件和協調條件。

(2)通過簡單剛柔耦合接觸系統的動力學仿真，驗證了動態粘合算法的計算精度和可靠性。

(3)復雜的接觸系統如鐵路車輛系統中受電弓與接觸網滑動接觸，輪軌滾動接觸等，使用動態粘合算法可以將這些系統分解成若干部件進行分別在不同的軟件或平臺上求解，通過數據傳輸和數據交換完成界面力和運動信息的傳遞，為接觸問題的分布式仿真提供了可行方案。

[1]Wang J Z，Ma Z D，Hulbert G M.A distributed mechanical system simulation platform based on a“gluing algorithm”[J].Journal of Computing and Information Science in Engineering，2005，5:71－76.

[2]Wang J Z.Development of advanced methodlogy for networkdistributed simulation[D].The University of Michigan，2005.

[3]Ryu G S.Integration of heterogeneous simulation models for network-distributed simulation [D]. The University of Michigan，2009.

[4]姚廷強，遲毅林，黃亞宇，等.剛柔耦合齒輪三維接觸動力學建模與振動分析[J].振動與沖擊，2009，28(2):167－171.

[5]周 寧;張衛華.基于直接積分法的弓網耦合系統動態性能仿真分析[J].中國鐵道科學，2008，129(16):71－76.

[6]王勖成.有限單元法[M].北京:清華大學出版社，2003.

[7]Kim H G.Interface element method(IEM)for a partitioned system with non-matching interfaces[J].Comput.Methods Appl.Mech.Engrg.，2002，191:3165 －3194.

[8]Kim H G.Interface element method:Treatment of nonmatching nodes at the ends of interfaces between partitioned domains[J].Comput.Methods Appl.Mech.Engrg.，2003，192:1841－1858.