脈沖子結構與有限元剛-彈混合連接的子結構方法

劉莉*,陳樹霖,周思達,陳昭岳

北京理工大學 宇航學院 飛行器動力學與控制教育部重點實驗室,北京 100081

隨著航天技術的發展以及對航天器性能要求的不斷提高,航天器結構變得日益復雜和龐大。在對現代復雜航天器進行動態響應分析以及結構優化時,采用傳統有限元方法對系統進行分析將引入眾多的自由度,消耗大量計算資源；面對系統歷經動響應分析的結構優化,傳統的有限元建模方法更是難以勝任。同時,航天器在研制過程中往往是不同部門、不同機構之間的合作,出于技術保護的考慮或者建模方式的不同,合作雙方無法直接共享模型。動態子結構方法在一定程度上能夠解決上述問題。動態子結構方法發展至今主要形成了兩類分支[1]:時域子結構方法和頻域子結構(FBS)方法。頻域子結構方法是在頻域內進行子結構綜合,1940年Biot第一次成功地應用該方法[2],之后又有不少學者不斷對頻域子結構方法進行發展完善[3-6]；時域子結構方法主要是指以Rayleigh-Ritz法為理論基礎的模態綜合,模態綜合從提出到現在已經得到了長足的發展[7-12],2010年Rixen[13]提出了一種新的時域子結構方法:基于脈沖響應函數(IRF)的子結構方法,簡稱脈沖子結構(IBS)方法。文獻[14]對該方法進行研究發展并將其應用于實際結構,結果表明該方法在處理小阻尼、瞬態沖擊問題時具有快速、高精度的優勢。之后不少學者在理論上對該方法進行擴展,文獻[15-16]對脈沖響應函數進行截斷處理并研究其對結果的影響,得到截斷的脈沖子結構方法；van der Valk和Rixen[17-18]研究了將基于脈沖響應函數的模型與有限元模型結合,但它們的連接方式僅僅考慮剛性連接,在航天器動響應分析中存在眾多剛-彈混合連接方式,比如月球探測器中心體與太陽翼之間的鉸鏈連接(仿真時可以將其視為平移自由度上的剛性連接與旋轉自由度上的彈性連接)；董威利等[19]研究了一種基于脈沖響應函數的彈性連接子結構綜合方法,該方法考慮了脈沖子結構間連接件的彈性特性,適合對連接元件進行優化設計,但方法未涉及脈沖子結構與有限元子結構間的連接,在進行系統結構優化時,需對局部結構進行修改,采用該方法每次優化迭代都將重新計算脈沖子結構的脈沖響應函數,而脈沖響應函數的多次重新生成將該方法高效的優勢大大削弱,使優化過程變得繁瑣耗時。

本文在上述研究的基礎上,發展了剛性、彈性以及剛-彈混合連接情形下脈沖子結構與有限元子結構的綜合方法,推導了相關的耦合計算格式,通過數值算例驗證了方法的正確性,并將該方法應用于月球探測器動響應分析,由于局部結構采用有限元建模,可以方便對其進行修改,該方法可被應用于大型系統動力學結構優化。

1 脈沖子結構方法及其降階形式

脈沖子結構方法于2010年由Rixen[13]第一次提出,它基于線性系統的疊加原理,基本思路是將子結構的脈沖響應函數與其所受全部外載荷在時域內進行卷積積分,通過界面相容條件將各子結構進行綜合得到總系統的位移響應,其數學本質是頻域子結構方法的傅里葉逆變換。

1.1 子結構動力學方程

將系統分成N個子結構,第s個子結構的運動方程為

式中:M(s)、C(s)和K(s)分別為第s個子結構的質量、阻尼和剛度矩陣；f(s)(t)和g(s)(t)分別為第s個子結構所受的外部激勵力列陣和界面處子結構受到的連接力列陣；¨u(s)(t)、˙u(s)(t)和u(s)(t)分別為第s個子結構的加速度、速度和位移列陣。

對于式(1)描述的線性系統,零初始條件下其動力學響應可以通過時域內Duhamel積分得到

式中:H(s)為子結構s的脈沖響應函數矩陣,其元素H(s)ij為子結構s的第j個自由度上受到單位脈沖函數δ(t)的激勵后,第i個自由度的位移響應。脈沖響應函數矩陣H(s)可通過直接積分算法(如Newmark法)得到也可以由試驗測得。

1.2 界面相容條件

整個系統的位移向量u可表示為

對于剛性連接的子結構有:連接界面處位移相等,連接力大小相等方向相反。通過布爾矩陣B將系統位移向量投影到子結構連接界面,位移相容性條件描述為

1.3 子結構運動綜合

綜合式(2)子結構動力學方程以及式(4)界面相容條件,同時利用式(5)將界面連接力通過拉格朗日乘子表達,得到全系統動力學方程為

該方程組即經典的脈沖子結構方法。

1.4 降階形式的脈沖子結構方法

子結構s為線性系統,由脈沖響應函數的物理實質可知,當第j個自由度上受激勵力(t)時,第i個自由度的位移響應(t)可由Duhamel積分得到

對于航天器,系統雖然復雜龐大,但設計人員關注的往往只是少部分關鍵位置處的響應輸出。同時,航天器的載荷邊界自由度相對較少,比如衛星與火箭的對接環、月球探測器與緩沖機構的連接點等。因此,對于第s個子結構只需保留界面自由度(t)、部分外載荷輸入自由度)(t)、以及除去界面自由度以外的關注位置處的自由度(t)。子結構s在)(t)自由度上進行激勵輸入,輸出自由度為(t),其表達式為

在降階后的子結構s中,重新構造布爾矩陣,使得

由式(7)～式(9),且子結構間通過剛性連接,得到降階后的系統動力學方程為

2 基于脈沖響應函數的子模型與有限元子模型的連接

由式(6)或式(10)可以得到多個脈沖子結構剛性連接的動響應結果,但在航天器結構中,有時不易獲得某些局部子結構的脈沖響應函數；在結構設計中,有時需要頻繁修改某局部結構,而通過直接積分求解脈沖響應函數過程繁瑣,這時將需修改的結構用有限元表示會更加方便。同時,航天器結構中,存在復雜的連接方式,僅僅通過剛性連接不能將其完全表示。接下來討論用不同連接方式將基于脈沖響應函數的子模型與有限元子模型進行連接綜合。

2.1 剛性連接

將整個系統分成兩類組成:基于脈沖響應函數的子模型和有限元子模型。系統動力學方程為

式中:上標s表示該變量屬于第s個脈沖子模型；上標r表示該變量屬于第r個有限元子模型。與第1節討論類似,f(t)和g(t)分別為外部激勵力列陣和界面處連接力列陣。

對于剛性連接,在界面處有

2.2 彈性連接

與剛性連接類似,容易得到彈性連接的界面相容條件為

式中:k為連接處的剛度矩陣。

由界面相容條件進而得到全系統的動力學方程為

接下來討論對該方程的離散求解。

對于式(16)的第1個式子采用梯形法則對其進行時間離散,有

記

式中:I為單位矩陣；k為彈性連接處的剛度矩陣。

將λn代入式(27)中的第1和第2個式子,可解得n時刻的位移,通過Newmark法進一步推進得到n時刻的加速度和速度響應。

2.3 剛-彈混合連接

連接界面可能存在剛-彈混合連接,比如月球探測器中心體與太陽翼的鉸接可以看成平動自由度上的剛性連接和旋轉自由度上的彈性連接。本節在之前討論的基礎上,進一步分析該連接模式下系統動力學方程的具體表達。

通過構造布爾矩陣B(sr)和B(sk)將脈沖子模型的位移un(s)投影到連接界面,有

其中:Irr和Ikk為單位矩陣,階數分別與連接處剛性連接和彈性連接的自由度數相同。

由式(33)容易看出,當界面全是彈性連接時,K退化成k,E退化成I；當界面全是剛性連接時,K退化成I,E退化成零矩陣。

將λn代入式(32)得到n時刻的位移,并通過Newmark法解得n時刻的加速度和速度響應。

3 數值算例

通過2組數值算例來驗證上述方法的可行性及正確高效性。第1組數值算例為簡單的彈簧質量塊系統,并將該系統劃分為兩個子結構,分別將兩子結構通過剛性、彈性、剛-彈混合連接進行組裝。在第2組數值算例中,該方法被應用到月球探測器動響應分析中,整個探測器被劃分為3部分,其中探測器中心體為基于脈沖響應函數的子結構,兩個太陽翼分別為兩個有限元子結構,太陽翼與探測器中心體通過剛-彈混合連接方式連接。兩組數值算例均與參考算例進行對比,參考算例通過Newmark積分(γ=0.5,β=0.25)得到。

圖1 彈簧質量塊系統示意圖(全彈性連接)Fig.1 Illustration of mass-spring system(substructures are coupled merely via springs)

3.1 彈簧質量塊系統

彈簧質量塊系統如圖1所示,整個系統共11個自由度,m1～m7為基于脈沖響應函數的子結構,m8～m11為有限元子結構,m5、m6、m7與m8可通過剛性、彈性、剛-彈混合3種方式連接,圖1中采用的是全彈性連接。系統各質量剛度參數如表1所示。

表1 彈簧質量塊系統參數Table 1 Parameters of mass-spring system

假設圖1中向右為正方向,在m9處施加矩形激勵力,如圖2所示。

圖2 激勵力示意圖Fig.2 Schematic of excitation force

以m10處的位移、速度和加速度結果為例,將數值算例結果與參考算例結果進行比對,參考算例通過Newmark積分(γ=0.5,β=0.25)得到。并采用相似模型誤差分析中常用的全局評價指標R2和RAAE[20]定量表示子結構方法的精度,R2越接近于1且RAAE越接近于0,表明預測結果的精度越高。

3.1.1 剛性連接

m5、m6、m7與m8之間均通過剛性連接,布爾矩陣為

m10處的位移對比如圖3所示。誤差評價指標結果如表2所示。

圖3 剛性連接時m10處的位移對比圖Fig.3 Comparison of displacement response of m10 when substructures are coupled merely by rigid joints

表2 剛性連接時m10響應的相對誤差Table 2 Relative error of response of m10 when substructures are coupled merely by rigid joints

3.1.2 彈性連接

m5、m6、m7與m8之間均通過彈性連接,布爾矩陣為

m10處的位移對比如圖4所示。誤差評價指標結果如表3所示。

圖4 彈性連接時m10處的位移對比圖Fig.4 Comparison of displacement response of m10 when substructures are coupled merely by springs

表3 彈性連接時m10響應的相對誤差Table 3 Relative error of response of m10 when substructures are coupled merely by springs

3.1.3 剛-彈混合連接

m5與m8之間通過剛性連接,而m6、m7與m8之間通過彈性連接,布爾矩陣為

m10處的位移對比如圖5所示。誤差評價指標結果如表4所示。

圖5 剛-彈混合連接時m10處的位移對比圖Fig.5 Comparison of displacement response of m10 when substructures are coupled by both springs and rigid joints

表4 剛-彈混合連接時m10響應的相對誤差Table 4 Relative error of response of m10 when substructures are coupled by both springs and rigid joints

由結果對比圖以及誤差評價表可看出,上述連接脈沖子結構模型與有限元模型的方法有效可行,并且位移、速度和加速度預測結果均保持很高精度。

3.2 月球探測器模型

本節使用的數值算例為腿式月球探測器模型,該模型如圖6所示。探測器由探測器中心體、著陸緩沖機構和有效載荷3部分組成[21]。足墊與緩沖腿、緩沖腿與探測器中心體通過萬向節連接,有效載荷按其相對位置以等效質量的形式分布于探測器中心體結構上。

圖6 月球探測器示意圖Fig.6 Configuration of lunar lander

探測器太陽翼采用復合材料鋪層。在結構設計及建模中,太陽翼的鋪層厚度、鋪層方向、鋪層分布等都需進行充分考慮。倘若采用優化算法對探測器太陽翼進行動力學結構優化,每次迭代必須將太陽翼與擁有眾多自由度的中心體組裝進行動響應分析,采用傳統的有限元方法進行響應預測將耗費大量計算資源。因此,很有必要引入子結構方法對原模型進行處理。

太陽翼屬于月球探測器的展開機構,它是通過類似扭簧的鉸鏈與中心體進行連接,連接形式是典型的剛-彈混合連接。

基于以上分析,將探測器進行子結構劃分:探測器中心體通過基于脈沖響應函數的子結構進行表征、兩個太陽翼均使用有限元模型建立。各太陽翼與探測器中心體均通過3個鉸鏈連接,每個鉸鏈在3個平動自由度上提供剛性約束,在轉軸方向上提供扭轉剛度。

太陽翼有限元模型如圖7所示,采用四節點完全積分殼單元,模型共包含3 024個自由度。

由于太陽翼質量遠小于探測器中心體質量(兩太陽翼質量不足探測器總質量的1%),其質量的改變對著陸過程中的緩沖作用力影響很小,在結構優化分析中,可利用探測器著陸沖擊數值試驗中測得的緩沖機構作用力作為激勵。激勵力作用于緩沖機構與中心體連接的12個連接點處。

考慮到子結構建模方式的不同,分別選取基于有限元的太陽翼角點和基于脈沖子結構的探測器中心體的頂板中心點作為比較對象(如圖6中所示的測點1和測點2)。

未減縮的探測器中心體有限元模型如圖8所示,共包含87 666個自由度,其中殼單元為四節點完全積分殼單元,涉及的梁單元采用二節點線性完全積分梁單元。按照降階形式的脈沖子結構的思路,保留激勵力輸入處36個自由度,連接界面處24個自由度,以及頂板中心點縱向的1個自由度。減縮后的探測器中心體,共60個輸入自由度,25個輸出自由度。

圖7 太陽翼有限元模型Fig.7 Finite element model of solar wing

圖8 探測器中心體有限元模型Fig.8 Finite element model of lander’s main body

在月球探測器以4 m/s的初速度、月球加速度g=1.63 m/s2的豎直著陸工況中測得緩沖機構對中心體的緩沖載荷。圖9給出了+Y方向主緩沖腿對中心體的縱向緩沖載荷隨著陸時間的變化曲線。

圖9 縱向緩沖載荷Fig.9 Longitudinal buffer load

使用Newmark積分方法(γ=0.5,β=0.25)對模型進行計算作為參考。+Y太陽翼角點處X方向的加速度響應結果如圖10所示。探測器中心體頂板處X方向的加速度響應結果如圖11所示。探測器兩測點X方向的動態響應相對誤差如表5所示。

圖10 太陽翼角點處的縱向加速度響應對比圖Fig.10 Comparison of longitudinal acceleration response of solar wing’s corner point

圖11 探測器中心體頂板測點處的縱向加速度響應對比圖Fig.11 Comparison of longitudinal acceleration response of lunar lander’s roof center

表5 子結構方法求得的探測器兩測點的動態響應相對誤差Table 5 Relative error of two test points’dynamic response solved by substructure method

在相同的計算條件下,參考算例計算耗時187.291 6 s,本方法計算耗時13.876 7 s,很大程度提高了計算效率；同時計算結果保持很高的精度,由于速度和加速度通過Newmark格式由位移遞推得到,加速度和速度的誤差略高于位移的誤差。

4 結 論

1)在經典脈沖子結構方法的基礎上,研究脈沖子結構和有限元子結構的剛性連接,提出了彈性以及剛-彈混合連接情形下的具體表達。通過彈簧質量塊算例,驗證了該方法的可行性和正確性。由于脈沖子結構保留了輸入輸出自由度上所有瞬態動力學特性,有限元子結構采用完整模型,該子結構綜合方法是一種保精度的動態子結構方法。

2)將本文所提方法應用于月球探測器結構動響應分析,用脈沖子結構方法建立的中心體模型與用有限元方法建立的太陽翼模型通過剛-彈混合方式連接。在保精度預測的前提下,很大程度上縮短了計算耗時。該方法適用于月球探測器著陸過程的動響應分析；太陽翼采用有限元建模,能夠方便對其進行參數修改,結合該方法的高效性,可以將其應用到探測器局部結構的動力學結構優化設計中。

[1] Schwarz H A.Gesammelte mathematische abhandlungen[M].Berlin:Springer-Verlag,1890:133-143.

[2] Biot M A.Coupled oscillations of aircraft engine-propeller systems[J].Journal of the Aeronautical Sciences,1940,7(9):376-382.

[3] Bishop R E D,Johnson D C.The mechanics of vibration[M].Cambridge:Cambridge University Press,1960:9-45.

[4] Rubin S.Mechanical immittance and transmission matrix concepts[J].The Journal of the Acoustical Society of America,1967,41(5):1171-1179.

[5] Jetmundsen B,Bielawa R L,Flannelly W G.Generalized frequency domain substructure synthesis[J].Journal of the American Helicopter Society,1988,33(1):55-64.

[6] Gordis J H.Structural synthesis in the frequency domain:A general formulation[J].Shock and Vibration,1994,1(5):461-471.

[7] Benfield A W,Hrude F R.Vibration analysis of structures by component mode[J].AIAA Journal,1971,9(7):1255-1261.

[8] Rubin S.Improved component-mode representation for structural dynamic analysis[J].AIAA Journal,1975,13(8):995-1006.

[9] Wang W L,Du Z R,Chen K Y.A short commentary for modal synthesis techniques and a novel improvement[J].Acta Aeronautica et Astronautica Sinica,1979(3):32-51(in Chinese).王文亮,杜作潤,陳康元.模態綜合技術短評和一種新的改進[J].航空學報,1979(3):32-51.

[10] Zheng Z C.Dynamic analysis of nonlinear systems by modal synthesis techniques[J].Applied Mathematics and Mechanics,1985,4(4):563-572(in Chinese).鄭兆昌.非線性系統動力分析的模態綜合技術[J].應用數學和力學,1985,4(4):563-572.

[11] Rixen D J.A dual craig-bampton method for dynamic substructuring[J].Journal of Computational and Applied Mathematics,2004,168(1-2):383-391.

[12] Dong W L,Liu L,Zhou S D.Model reduction of soft landing dynamics for lunar lander with local nonlinearities[J].Acta Aeronautica et Astronautica Sinica,2014,35(5):1319-1328(in Chinese).董威利,劉莉,周思達.含局部非線性的月球探測器軟著陸動力學降階分析[J].航空學報,2014,35(5):1319-1328.

[13] Rixen D J.A substructuring technique based on measured and computed impulse response functions of components[C]//The 24th International Conference on Noise and Vibration Engineering(ISMA 2010).Leuven,Belgium:Katholieke Universiteit Leuven,2010:1939-1954.

[14] Dong W L,Liu L,Zhou S D.Acceleration response prediction for lunar lander using time-domain substructure methods[J].Acta Aeronautica et Astronautica Sinica,2014,35(3):848-856(in Chinese).董威利,劉莉,周思達.月球探測器加速度響應預測的時域子結構方法[J].航空學報,2014,35(3):848-856.

[15] Rixen D J,Haghighat N.Truncating the impulse responses of substructures to speed up the impulse-based substructuring[C]//30th IMAC.New York:Springer,2012:137-148.

[16] van der Seijs M V,Rixen D J.Efficient impulse based substructuring using truncated impulse response functions and mode superposition[C]//International Conference on Noise and Vibration Engineering 2012(ISMA 2012).Leuven,Belgium:Katholieke Universiteit Leuven,2012:3487-3499.

[17] van der Valk P L C,Rixen D J.An effective method for assembling impulse response functions to linear and nonlinear finite element models[C]//30th IMAC.New York:Springer,2012:123-135.

[18] van der Valk P L C,Rixen D J.An impulse based substructuring method for coupling impulse response functions and finite element models[J].Computer Methods in Applied Mechanics and Engineering,2014,275:113-137.

[19] Dong W L,Liu L,Zhou S D,et al.Substructure synthesis in time domain with rigid-elastic hybrid joints[J].AIAA Journal,2015,53(2):504-509.

[20] Jin R,Chen W,Simpson T W.Comparative studies of metamodelling techniques under multiple modelling criteria[J].Structural and Multidisciplinary Optimization,2001,23(1):1-13.

[21] Otto O R,Laurenson R M,Melliere R A,et al.Analyses and limited evaluation of payload and legged landing system structures for the survivable soft landing of instrument payloads,NASA CR-111919[R].Washington,D.C.:NASA,1971.