針對局部非線性問題的混合坐標模態綜合法

王陶，何歡，2，＊，陳國平，2

1．南京航空航天大學 機械結構力學及控制國家重點實驗室，南京 210016 2．南京航空航天大學 振動工程研究所，南京 210016

針對局部非線性問題的混合坐標模態綜合法

王陶1，何歡1，2，＊，陳國平1，2

1．南京航空航天大學 機械結構力學及控制國家重點實驗室，南京 210016 2．南京航空航天大學 振動工程研究所，南京 210016

非線性動力系統模型的計算效率問題是結構動力學領域中的重要研究課題。提出了一種針對局部非線性問題的混合坐標自由界面子結構模態綜合方法。根據局部非線性系統的特點，將結構按照線性部分與非線性部分進行分割。線性部分子結構可以通過模態坐標轉換到模態空間。在對線性部分進行減縮的過程中考慮了剩余柔度的影響，并通過構造一組與低階模態關于系統矩陣加權正交的向量組，解決了子結構含有剛體模態時剩余柔度矩陣無法計算的問題。非線性部分子結構則保留原有的物理坐標。通過界面協調關系，采用自由界面方法得到系統混合坐標綜合方程。最后，通過數值算例驗證了所提出方法的有效性。

局部非線性；混合坐標；自由界面；模態綜合；協調關系

隨著工程結構的日益復雜化，研究者越來越重視結構的非線性特征對結構動力學行為的影響。其中有一類結構系統在變形過程中只有局部呈現出明顯的非線性特征，而其余部分仍保持線性，這就是所謂的局部非線性問題。例如，飛行器著陸沖擊過程［1］就是一個典型的局部非線性問題。對于這類局部非線性問題，如果采用有限元法（FEM）對其直接進行建模，所得到的動力學方程為非線性動力學方程，當對其進行動力學分析時，其計算量將會數倍于線性模型。于是，在保證計算精度的基礎上如何有效降低此類問題的計算規模就成為許多研究者關注的一個重點。其中，子結構模態綜合法就是一種減縮結構自由度數目的經典方法。

子結構模態綜合法最早由 Hurty［2］于1965年提出，文中主要討論了無阻尼線性結構的模型減縮方法。Goldman［3］和 Macneal［4］針對比例阻尼結構提出了相應的改進。隨后，很多學者在提高計算精度、處理復雜界面條件等方面也進行了大量的研究工作［5－7］。在不斷的改進與發展之后，逐漸形成了目前針對線性結構的較為完備的理論體系。近年來，越來越多的研究者開始將子結構模態綜合法應用于更為復雜的結構動力學問題中去，例如在一般黏性阻尼結構振動分析［8］、大型復雜結構模型修正［9］、結構不 確定性 分析［10－11］等問題中模態綜合法都表現出了很好的應用前景。

對于局部非線性問題來說，非線性主要存在某一小范圍區域內，其他大部分結構可以認為仍然符合線彈性假設。利用這個特點，很多學者也開始把子結構模態綜合法的思想引入局部非線性結構的動力學分析中，通過對模型進行減縮來提高其計算效率。Clough－Wilson［12］最先提出將子結構模態綜合法應用于大型局部非線性結構的動力學分析當中。郝淑英等［13］采用子結構模態綜合法建立了非線性連接元組成的連接子結構的減縮模型，并通過漸進法求解了系統的一次近似定常解。文獻［14－16］采用子結構模態綜合法對轉子系統的局部非線性問題進行了研究。Verros和Natsiavas［17］用自由界面模態綜合法對機車的非線性動力學問題進行了相關研究。Kawamura等［18］提出了一種針對多自由度非線性系統的子結構模態綜合法，通過計算結果可以看出在子結構沒有剛體模態時，其方法具有較高的計算精度。Saito等［19］采用 Craig－Bampton方法對裂紋結構的非線性受迫振動進行了動響應分析。Praveen Krishna和Padmanabhan［20］采用不同模型減縮方法對局部非線性系統進行了分析，并對這些減縮方法的有效性進行了討論。上述研究在采用自由界面子結構模態綜合法對局部非線性問題進行分析時，沒有考慮剩余柔度的影響，且在子結構存在剛體模態時的計算精度較差。

針對上述問題本文提出了一種混合坐標自由界面子結構模態綜合法。為了提高計算精度，考慮了剩余柔度的影響，并解決了含剛體模態時的系統剩余柔度矩陣無法直接求解的問題。在綜合過程中對線性部分進行模態減縮，非線性部分仍保持物理坐標下的形式，所以系統依然保留了原有的非線性特性。通過界面協調條件將轉換到模態坐標下的線性子結構與物理坐標下的非線性子結構組裝得到由模態坐標與物理坐標混合表示出的綜合方程。由于綜合模型規模的減小使得其計算效率可以得到很大的提高。

1 線性子結構模型減縮

根據局部非線性系統的特點，將系統按照線性結構與非線性結構進行劃分，將線性子結構記為子結構a，則其結構動力學方程可表示為

式中：Ma、Ca和Ka∈Rn×n分別為子結構a的質量矩陣、阻尼矩陣和剛度矩陣；ua和fa∈Rn×1分別為廣義位移向量和載荷向量。通過特征值分析，容易計算出線性子結構較低階的l階模態Φla。需要指出的是，由于考慮自由界面形式，所以這里Φla可能包含有剛體模態。

1．1 構造加權正交向量集

假設存在矩陣φa∈Rn×（n－l）滿足

使得

先給定線性無關的n－l個向量

結合式（2）和式（3）可解得

由此可得

1．2 剩余柔度矩陣與界面協調方程

定義Φa＝［Φlaφa］，則子結構的物理坐標可以表示為

利用Φa可將式（1）轉換為

將式（9）展開，其高階模態坐標所對應的方程可以表示為

將式（13）進行Laplace變換后，可得式中：s為復變量；F珚a為fa在s域中的表示形式。

將式（14）進行泰勒展開，并取其一階近似，再進行反Laplace變換后，可得

這實際上相當于保留了剛度項的影響。

將式（15）代入式（8），展開后可得

式中：

式中：

其中：Ga∈Rn×n為剩余柔度矩陣。由于φa為構造出來的等效高階模態，φa中并未包含剛體模態。所以式（9）與式（12）中的僅與高階模態有關，為可逆矩陣。因此，當子結構存在剛體模態時，仍可以通過式（17）來求解剩余柔度矩陣。

假設子結構的交界面上有m個自由度，將廣義物理坐標進行分塊處理，按照內部坐標與界面坐標寫為

式中：uaj為子結構界面自由度所對應的位移；uai為子結構內部自由度所對應的位移；fai為結構內部自由度上所受到的外力；faj為界面力為界面位移所對應的界面模態；為內部位移所對應的內部模態。

將式（18）展開，則界面坐標可以表示為

2 非線性子結構動力學方程

將非線性子結構記為子結構b，其動力學方程可表示為式中：Mb、Cb和Kb分別為子結構b線性部分的質量矩陣、阻尼矩陣和剛度矩陣；g（ub， ub）為其非線性項。

類似地，將其按照內部自由度與界面自由度進行分塊，可以寫為

整理式（21），子結構b的界面位移可以表示為

3 局部非線性系統的混合坐標綜合方程

僅考慮簡單的界面情況，即界面外載荷為零，界面上不添加集中質量和彈簧，所以界面位移的連續性條件與載荷協調關系可寫為

式中：

將式（25）代入式（9）與式（21），可以得到結構的綜合方程為

式中：

定義P＝ （Φlaj）TGabK－1bjj，則式（28）～式（30）中的相關參數表達式為

4 數值算例

4．1 算例1

如圖1所示為一個考慮阻尼的十自由度彈簧質量塊系統（阻尼器沒有出現在圖中，但認為阻尼器平行于彈簧存在于系統中），其中最右側與固定端連接的彈簧具有非線性特性，這里非線性力可以表示為fnonlinear＝k（x＋εx3）。將系統在第10個質量塊處分割為兩個子結構，記左側結構為子結構1，其僅具有線性特性；記右側結構為子結構2，其具有非線性特性。系統的相關參數為

式中：γj為Reyleigh阻尼常數；ε為小參數。

圖1 十自由度局部非線性系統Fig．1 A ten－degree－of－freedom system with localised nonlinearities

分別提取子結構1的2、4、6和8階低階模態進行綜合，對非線性系統所對應的派生系統進行特征值分析。表1比較了綜合模型與原模型計算得到的固有頻率，其中r為保留的低階模態數。從表1中可以看出，隨著提取的低階模態數目的增加，通過本文方法計算所得到的固有頻率的精度也逐漸提高。

假設一個周期外力作用于m9，采用諧波平衡法對非線性系統進行頻響分析。其中周期力可表示為：f9＝F9cos（ωt＋θ），F9＝500N，其中ω為激勵頻率，θ為激勵與響應之間的相位差。若只考慮一階諧波，則根據諧波平衡法此時系統的穩態響應可以表示為：xj＝Ajcos（ωt）（j＝1，2，…，l＋1），其中Aj為響應幅值。將周期力及響應的表達式代入綜合方程，通過比較等式兩邊的一次諧波項系數，就可以得到相應的幅頻特性曲線表達式。

如圖2所示，為系統第一階固有頻率所對應的幅頻特性曲線。從圖中可以看出，系統的幅頻特性曲線有著明顯的跳躍現象，且隨著提取的低階模態數目的增加，通過模態綜合法計算所得的幅頻特性曲線也逐漸趨近于原始結構的計算結果。

圖2 m10的幅頻特性曲線比較Fig．2 Comparison of frequency－amplitude curves of m10

4．2 算例2

考慮如圖3所示的具有局部非線性的飛機著陸問題的結構系統。將其分為機身部分與起落架部分。其機身部分為線性區域，作為子結構a，通過用空間梁單元對其進行建模，共有69個梁單元，414個自由度。機身部分的結構阻尼假設為：C＝αM＋βK，α與β為Rayleigh阻尼常數。機身結構的總質量為9 195．75kg。機身模型的其他參數如表2所示。

對于具有非線性特性的起落架系統（如圖3中紅色圓圈所示為其安裝位置）通過引入非線性彈簧阻尼單元來進行模擬。此部分作為子結構b。起落架內部的緩沖力可以表示為fnonlinear（u， u）＝kbu＋γkbu3＋c1u。其中u為沿起落架軸線方向的位移。采用集中質量單元模擬起落架系統的質量，每個起落架的等效質量記為m1。起落架模型的參數如表3所示。

以圖3所示結構在平衡位置處的微振動為基準。在此基礎上提取通過模態分析得到的子結構a的前10、20、30和40階低階模態（包含前6階剛體模態）進行模態綜合。表4分別給出了完全有限元模型與模態綜合模型的計算結果。

圖3 飛機總體結構模型Fig．3 Model of fuselage structure

表2 機身模型參數Table 2 Parameters of fuselage model

表3 起落架模型參數Table 3 Parameters of landing gears model

表4 模態綜合模型與完全有限元模型的固有頻率比較Table 4 Comparison of natural frequencies between synthesis models and full FE model

圖4 固有頻率誤差對比Fig．4 Comparison of errors for natural frequencies

圖4所示為保留前20、30和40階低階模態所計算得到的前20階固有頻率的誤差對比。從圖4中可以看出，通過本文方法所計算得到的固有頻率有較高的精度，當提取低階模態數目達到40階時，計算的誤差都在1%以內。考慮機身在z方向上初始速度為－5m／s時的動力學時域響應分析。分別采用子結構綜合模型和完全有限元模型進行計算，仿真分析時間為3s。采用本文提出的模態綜合法，分別保留線性子結構的前10、20、30和40階低階模態用于模態綜合。

圖5給出了圖3結構上Point 1、Point 2與Point 4沿z方向上的位移響應曲線比較。從圖5可以看出，位移響應與完全有限元模型計算結果吻合得很好。隨著采用的低階模態數目的增加，運用模態綜合所計算出的位移響應都與完全有限元模型給出的計算結果趨于一致。

圖5 Points 1，2，4沿z方向的位移響應Fig．5 Displacement response along z direction at Points 1，2，4

定義CPU time ratio為Ts／Tf，其中Ts為綜合模型計算所需CPU計算時間，Tf為完全有限元模型計算所需CPU計算時間。從表5中可以看出，綜合模型的計算效率有了非常大的提高。

表5 綜合模型與完全有限元模型計算時間比較Table 5 Comparison of CPU time consumed between synthesis models and full FE model

5 結 論

1）在綜合過程中考慮了剩余柔度的影響，通過構造一組與低階模態關于系統矩陣加權正交的向量組，解決了子結構具有剛體位移時剩余柔度矩陣無法求解的問題。采用自由界面綜合方法得到了混合坐標下的模態綜合方程。

2）利用諧波平衡法對一多自由度彈簧質量塊系統進行了頻響分析。計算結果表明本文方法可以有效地處理局部非線性系統頻域下的振動特性問題，這為進一步研究諸如分岔等非線性行為提供了一種實用方法。

3）對飛機著陸問題進行了時域響應分析。數值分析結果表明，利用本文方法得到的減縮模型能夠高效準確地求解局部非線性系統的動響應問題。合理保留一定數量的低階模態后，綜合模型不僅具有良好的計算精度，更重要的是從本質上改善了局部非線性系統動響應問題的計算效率。

［1］ 董威利，劉莉，周思達．含局部非線性的月球探測器軟著陸動力學模型降階分析［J］．航空學報，2014，35（5）：1319－1328．DONG W L，LIU L，ZHOU S D．Model reduction analysis of soft landing dynamics for lunar lander with local nonlinearities［J］．Acta Aeronautica et Astronautica Sinica，2014，35（5）：1319－1328（in Chinese）．

［2］ HURTY W C．Dynamic analysis of structural systems using component modes［J］．AIAA Journal，1965，3（4）：678－785．

［3］ GOLDMAN R L．Vibration analysis by dynamic partitioning［J］．AIAA Journal，1969，7（6）：1152－1154．

［4］ MACNEAL R H．A hybrid method of component mode synthesis［J］．Computers and Structures，1971，1（4）：581－601．

［5］ RUBIN S．Improved component－mode representation for structural dynamic analysis［J］．AIAA Journal，1975，13（8）：995－1006．

［6］ QIU J B，YING Z G，WILLIAMS F W．Exact modal synthesis techniques using residual constraint modes［J］．International Journal for Numerical Methods in Engineering，1997，40（13）：2475－2492．

［7］ BENFIELD W A，HRUDE R F．Vibration analysis of structures by component mode substitution［J］．AIAA Journal，1971，9（7）：1255－1261．

［8］ HE H，WANG T，CHEN G P，et al．A real decoupled method and free interface component mode synthesis methods for generally damped systems［J］．Journal of Sound and Vibration，2014，333（2）：584－603．

［9］ PAPADIMITRIOU C，PAPADIOTI D C．Component mode synthesis techniques for finite element model updating［J］．Computers and Structures，2013，126（1）：15－28．

［10］ CHENTOUF S A，BOUHADDI N，LAITEM C S．Robustness analysis by aprobabilistic approach for propagation of uncertainties in a component mode synthesis context［J］．Mechanical Systems and Signal Processing，2011，25（7）：2426－2443．

［11］ MENCIK J M．Model reduction and perturbation analysis of wave finite element formulations for computing the forced response of coupled elastic systems involving junctions with uncertain eigenfrequencies［J］．Computer Methods in Applied Mechanics and Engineering，2011，200（45－46）：3051－3065．

［12］ CLOUGH R W，WILSON E L．Dynamic analysis of large structural systems with local nonlinearity［J］．Computer Methods in Applied Mechanics and Engineering．1979，17（18）：107－129．

［13］ 郝淑英，陳予恕，張琪昌，等．連結子結構在非線性動力學分析中的應用［J］．天津大學學報，2001，34（3）：295－299．HAO S Y，CHEN Y S，ZHANG Q C，et al．Application of link substructure to nonlinear dynamic system analysis［J］．Journal of Tianjin University，2001，34（3）：295－299（in Chinese）．

［14］ 華軍，許慶余，張家忠．擠壓油膜阻尼器－滑動軸承－轉子系統非線性動力特性的數值分析及實驗研究［J］．航空學報，2001，22（1）：42－45．HUA J，XU Q Y，ZHANG J Z．Numerical and experimental study on nonlinear dynamic behavior of the fluid film bearing－rotor system with squeeze film damper［J］．Acta Aeronautica et Astronautica Sinica，2001，22（1）：42－45（in Chinese）．

［15］ IWATSUBO T，SHIMBO K，KAWAMURA S．Nonlinear vibration analysis of a rotor system using component mode synthesis method［J］．Archive of Applied Mechanics，2003，72（11－12）：843－855．

［16］ 呂延軍，虞烈，劉恒．非線性轉子－軸承系統的動力學特性及穩定性［J］．機械強度，2004，26（3）：242－246．LV Y J，YU L，LIU H．Dynamic characterstics and stability of nonlinear rotor－bearing system［J］．Journal of Mechanical Strength，2004，26（3）：242－246 （in Chinese）．

［17］ VERROS G，NATSIAVAS S．Ride dynamics of nonlinear vehicle models using component mode synthesis［J］．Journal of Vibration and Acoustics，2002，124（3）：427－434．

［18］ KAWAMURA S，NAITO T，ZAHID H M，et al．Analysis of nonlinear steady state vibration of a multi－degreeof－freedom system using component mode synthesis method［J］．Applied Acoustics，2008，69（7）：624－633．

［19］ SAITO A，CASTANIER M，PIERRE C，et al．Efficient nonlinear vibration analysis of the forced response of rotating cracked blades［J］．Journal of Computational ＆ Nonlinear Dynamics，2006，4（1）：53－63．

［20］ PRAVEEN KRISHNA I R，PADMANABHAN C．Improved reduced order solution techniques for nonlinear systems with localized nonlinearities［J］．Nonlinear Dynamics，2011，63（4）：561－586．

Component mode synthesis method based on hybrid coordinates for structure with localised nonlinearities

WANG Tao1，HE Huan1，2，＊ ，CHEN Guoping1，2
1．State Key Laboratory of Mechanics and Control of Mechanical Structures，Nanjing University of Aeronautics and Astronautics，Nanjing 210016，China 2．Institute of Vibration Engineering Research，Nanjing University of Aeronautics and Astronautics，Nanjing 210016，China

The calculation efficiency of nonlinear dynamic system has become increasingly important in the structural dynamics field．A hybrid coordinates component mode synthesis method is proposed in this paper for the structure with localised nonlinearities．According to its feature，generally，the system is divided into the linear component and nonlinear component．The equations of the linear component can be transformed into the modal coordinates by its linear vibration modes．In order to improve the accuracy，the residual flexibility attachment matrix of the system is introduced．And by constructing the weighted－orthogonal vector sets which have weighted－orthogonal relationship with the lower retained modes，the residual flexibility attachment matrix is obtained without using inverse of the stiffness matrix．The equations of the nonlinear component are kept as their original form．The synthesis equations which are expressed by hybrid coordinates are derived in terms of compatibility conditions at the interface．Finally，applications of the proposed methods to the numerical examples demonstrate that the present method is computationally effective．

localised nonlinearities；hybrid coordinates；free interface；component mode synthesis method；compatibilityconditions

2015－09－09；Revised：2015－12－08；Accepted：2015－12－28；Published online：2016－01－20 13：47

URL：www．cnki．net／kcms／detail／11．1929．V．20160120．1347．002．html

s：National Natural Science Foundation of China（11472132）；the Fundamental Research Funds for Central Universities（NS2014002）；Priority Academic Program Development（PAPD）of Jiangsu Higher Education Institutions

V215；O322；TB123

A

1000－6893（2016）09－2757－09

10．7527／S1000－6893．2015．0360

2015－09－09；退修日期：2015－12－08；錄用日期：2015－12－28；網絡出版時間：2016－01－20 13：47

www．cnki．net／kcms／detail／11．1929．V．20160120．1347．002．html

國家自然科學基金（11472132）；中央高校基本科研業務費專項資金（NS2014002）；江蘇高校優勢學科建設工程

＊通訊作者．Tel．：025－87138242 E－mail：hehuan＠nuaa．edu．cn

王陶，何歡，陳國平．針對局部非線性問題的混合坐標模態綜合法［J］．航空學報，2016，37（9）：27572－765．WANG T，HE H，CHEN G P．Component mode synthesis method based on hybrid coordinates for structure with localised nonlinearities［J］．Acta Aeronautica et Astronautica Sinica，2016，37（9）：27572－765．

王陶 男，博士。主要研究方向：復雜結構動力學。Tel．：025－84892197

E－mail：wangtao813619＠163．com

何歡 男，博士，副教授，碩士生導師。主要研究方向：復雜結構動力學、飛行器回收系統動力學。

Tel．：025－87138242

E－mail：hehuan＠nuaa．edu．cn

陳國平 男，博士，教授，博士生導師。主要研究方向：復雜結構動力學與控制。

Tel．：025－84892142

E－mail：gpchen＠nuaa．edu．cn

＊Corresponding author．Tel．：025－87138242 E－mail：hehuan＠nuaa．edu．cn