一種多體系統沖擊響應譜分析法

賀少華，吳新躍

(海軍工程大學 船舶與動力學院，武漢 430033)

對于結構動力學，當然也包括運動學，響應特性的獲得一般分為三個步驟，首先建立實際對象的數學模型由于實際對象理所當然是一個連續系統，所以對于簡單結構對象，我們可以直接以偏微分方程的形式建立連續系統數學方程，但是我們對這種方程的求解能力十分有限，一般僅限于弦、桿、軸、梁等簡單結構。當對象結構變得復雜時，我們必須采用另外一種建模方法，那就是將連續系統離散化，即建立連續系統的“近似”—離散系統。目前建立離散系統的方法有很多，如我們熟悉的集中質量法、有限元方法、多體系統法等。建立完實際對象的抽象數學模型，之后就是要建立模型在特定物理環境中的動力學方程，目前動力學方程的建立主要基于兩大力學理論——矢量力學和分析力學，基于這兩大理論，出現了一大批風格各異各有所長的動力學方程建立方法，如牛頓－歐拉方程法、拉格朗日方程法、變分方法、凱恩方法、旋量方法、傳遞矩陣方法等，其中的很多方法是隨著計算機技術的發展而建立的，其共同目標是要實現一種高度程式化，適應編制計算程序的動力學方程建立方法。最后，得到物理模型的動力學方程后，就是對這些方程的求解，從而得到動力學響應特性。對于線性系統方程而言，我們可以采用解析方法、模態疊加法等，對于非線性方程而言，我們則可以采用近似解析方法和數值積分方法，近似解析方法主要針對弱非線性系統的穩態響應的求解，如攝動法、諧波平衡法、多尺度方法、平均法等，而數值積分方法如顯式算法中的中心差分法、Runge_Kutta法和隱式算法中的 Wilson-θ法、Houbolt法、Newmark法等，它們則適合強非線性問題的求解［1－4］。當然，上面的這些方法有些并沒有嚴格的界限，有些方法的本質是一致的。另外，還有許多方法是在這些方法基礎之上或者聯合幾種方法建立來的，在這就不一一列舉。

艦用機械設備沖擊響應仿真建模與計算是隨著計算機軟硬件水平和力學、計算數學等學科知識的發展而發展起來的，幾十年來，已出現了數十種方法，均體現了當時相應各學科的發展水平，總的來說，這些方法均有特定適應的對象，各有優缺點，針對不同的分析對象必須采用不同的計算分析方法。顯然，評價一個力學模型優劣的重要標準應該是該模型是否能夠可靠和高速處理各種動力學現象，通常解的精確和計算所要付出的代價是一對矛盾。對于艦用設備特別是復雜動力機械設備，通常采用有限元法、有限元子結構法、多體系統法、集總質量法進行離散化建模。而沖擊動力學方程的建立則主要采用經典矢量力學的方法，方程響應的求解則主要采用模態疊加法、數值積分法。針對線性系統，主要采用基于模態疊加原理的譜分析法［5－7］，如著名的 DDAM 法就是其中的一種;針對非線性系統，則主要出現了增量模態疊加法、偽力法和幾種數值積分方法［8，9］。

隨著有限元和計算機技術水平的飛速發展，目前借助商業化的有限元軟件進行艦用機械設備的沖擊響應計算分析成為了國內的主流。但是，必須指出的是，絕大部分有限元商業軟件的知識產權屬于國外，過分依賴這些商業軟件將使我國本已滯后的艦船及設備抗沖研究特別是理論研究更加落后于國外。實際上，在上面文字中已經提到，除有限元方法外，另外一種建模和計算方法即多體系統法在艦用復雜機械沖擊動力學研究上更具應用潛質，因為該方法相對有限元法具有更加節省計算資源和易于理論創新的優勢。近幾年來，國內相關學者開始將多體系統動力學方法運用到艦船及設備抗沖研究上面。我們知道，目前的沖擊響應譜分析方法一般采用的還是有限元建模，而本文的目的就是要建立一種基于多體系統建模方法的沖擊響應譜分析方法。

1 方法理論基礎

我們知道，沖擊響應譜分析方法的基礎是模態疊加理論，系統特征矢量(特征函數)的正交性是系統動力響應可以用有限幾階模態就可快速收斂到所需工程精度要求的前提條件，但是，由于多體系統中剛體與彈性體之間的動力耦合作用，使得多剛柔系統的特征值問題非自共軛，振型函數不具有通常意義下的正交性，難以用模態方法精確分析系統的動力響應。為解決多剛柔系統特征矢量正交性，文獻［10，11］引入了多體系統增廣特征矢量和增廣算子的概念，通過建立滿足對稱性的增廣公式，證明了增廣特征矢量自動滿足正交性，為實現多剛柔系統動力響應的精確分析創造了條件，也為本文的多體系統沖擊響應譜分析法的創立提供了可能性。增廣特征矢量同時含有體元件線位移和角位移的模態坐標，但又不是兩者的簡單結合，同時含有離散變量和連續變量，它的變量個數等于系統體動力學方程的個數，這些都是增廣特征矢量不同于一般特征矢量的特點。增廣特征矢量的構造過程簡單而程式化，只需按序列寫出描述體元件運動狀態的位移參量即可得到增廣特征矢量。

下面我們以一個基礎受沖擊簡單多剛柔系統為例來闡述本文提出的多體系統沖擊響應譜分析法的一般原理。

如圖1所示平面振動剛柔系統，由彈性系數為k的彈簧1、質量為m的集中質量(剛體)2和長度為l、線質量密度為、抗彎剛度為EI等截面Euler-Bernoulli梁3組成。集中質量2受光滑導槽的約束做往復運動，梁3左端與剛體2固接于點Q，右端自由。初始時刻梁3未變形，基礎的y向加速度沖擊激勵時間歷程為ü，無外力和阻尼作用，求系統任意點沖擊響應時間歷程。

圖1 由彈簧、集中質量和梁組成的剛柔混合系統Fig.1 A multi-rigid-f lexible-body system composed of spring，centralized mass and beam

有限元法是結構分析中應用最廣泛和強有力的工具，但對復雜結構使用大量節點導致必須處理和管理非常大的矩陣，計算工作量大，大剛度梯度系統往往伴隨計算“病態”，導致計算困難，而多體系統結合傳遞矩陣法就能較好地克服這一問題。

以彈簧1與地面(基礎)0的聯結點為輸入端，梁的自由端 4 為輸出端，定義狀態矢量 Z0，Z1，2，Z2，3，Zx(0 ＜ x≤l)為 Z=［Y，Θz，Mz，Qy］T，其中 Y，Θz，Mz，Qy分別為沿y軸向的位移，繞z軸的角位移，z向力矩和y向剪切力。各個元件的傳遞方程為:

各個元件的傳遞矩陣為:

系統的總傳遞方程為:

將邊界條件:

代入總傳遞方程，得特征方程:

數值求解上式可得系統固有頻率 ωk(k=1，2，…)。對每一階ωk，求解總傳遞方程可得狀態矢量Z0，再由各個元件之間的傳遞方程可得系統任意點的狀態矢量和系統振型。

建立多剛柔系統的運動微分方程:用y3(x，t)(0＜x≤l)表示梁上P點的位移，y2，1(t)表示集中質量2的位移，則系統的運動微分方程為:

將上述方程組轉化為相對基礎坐標系中的方程組:

從上兩式可以發現，基礎加速度沖擊激勵已很巧妙地引入到了系統動力學方程中。

將上面動力學方程組寫成體動力學方程的形式:

式中:

算子D3的物理意義是:當D3作用于系統中某聯接點的位移時，表示體元件在該聯接點處受到的除阻尼力外的所有內力，“I”表示輸入端受到的內力，“O”表示輸出端受到的內力。

定義模態參與因子pj為:

對于零初始條件，其模態位移的解(模態系數或廣義坐標)為:

沖擊速度譜由如下卷積定義［5］:

于是，最大可能的模態位移的解(模態系數或廣義坐標)為:

因此，模態j的最大可能位移為:

求得沖擊激勵下每個振動模態最大響應后，需要對模態進行組合。在時域內，各模態貢獻的代數和導致諸質量點的物理位移，可考慮相位和符號。然而，可使用各種方法來組合最大模態撓度，通常可使用下列三種方法之一:

① 絕對值求和(ABS):

② 平方和的均方根(SRSS):

③ 美國海軍研究實驗室求和(NRL):

2 工程應用實例

以圖1所示系統為例，取 K=615.86 N/m，m=15.6kg，l=5 m=6.24kg/m，EI=213.33 N·m2，用三種方法計算得出的固有頻率如表1所示。顯然，用本文中的多體系統結合傳遞矩陣方法得到的系統的振動特性是精確解，有限元方法的解是近似解。

表1 固有頻率計算結果/(rad·s－1)Tab.1 Natural frequencies/(rad·s－1)

系統所受沖擊速度譜如圖2所示，用上述參數，三種方法(解析方法、有限元法和本文方法)計算出的點4的位移響應時間歷程如圖3所示，取前7階模態響應進行NRL方法合成。其他參數不變，取K=10000 N/m，計算結果如圖4所示。可見，無論哪一種情況，本文方法得出的解(圖3、圖4中“△”表示)始終與解析解基本重合(圖3、圖4中粗實線)，而有限元方法在第一種情況時與解析解重合較好(圖3細實線)，當m減小或是k增大時，有限元法的誤差變大(圖4細實線)。必須說明的是，這里所謂的解析解指的是按一般振動力學方法建立振動方程，然后根據響應譜分析法(模態疊加法)得到系統的響應，疊加的時候同樣選用NRL合成方法。對于該工程實例來說，由于結構簡單，所以可以計算得到解析解(當然，這里的“解析”也是一種近似解，因為截取的模態數有限)，在計算中，我們發現，即使是如此簡單的結構，解析解的表達式還是很復雜的，而本文方法只需程式化地進行元件的傳遞矩陣拼裝，即可方便地計算出系統的固有振動特性，且固有振動特性對應于精確解析解，所以求解的動力響應精度高，這就出現了在該工程實例中解析解和本文方法計算得到的解幾乎重合的情況。我們知道，對于一般的復雜系統，找到解析解是不可能的，這時候，本文方法的優越性更能得以體現。

圖2 沖擊輸入譜Fig.2 Shock input spectrum

圖3 梁自由端的位移時間歷程Fig.3 Time-varying displacement of free end of beam

圖4 k=10000 N/m時梁自由端位移時間歷程Fig.4 Time-varying displacement of free end of beam when k=10000 N/m

3 結論

綜合全文內容，我們可以發現，結合傳遞矩陣法的多體系統沖擊響應譜分析法具有以下特點:

(1)涉及的矩陣階次不取決于系統的自由度數，所以計算量比通常方法要小得多，如果對本文實例一類的鏈式系統，涉及的矩陣階次僅取決于元件的矩陣階次，這對工程計算至關重要;

(2)對于連續系統，其運動微分方程為偏微分方程，但本文方法只需求解常微分方程即可得到系統的沖擊響應;

(3)本文方法可較準確地獲得多體系統的真實振型和沖擊響應時間歷程。

［1］劉延柱，洪嘉振，楊海興.多剛體系統動力學［M］.上海:高等教育出版社，1989.

［2］Kim S J，Cho J Y，Kim W D.From the trapezoidal rule to higher-order accurate and unconditionally stable time-integration method for structural dynamics［J］.Computer Methods in Applied Mechanics and Engineering，1997，149:73 －88.

［3］Fung T C.Unconditionally stable higher-order accurate collocation time-step integration algorithms for first-orderequations［J］.Computer Methods in Applied Mechanics and Engineering，2000，190:1651 －1662.

［4］Lee A S，Kim B O，Kim Y C.A finite element transient response analysis method of a rotor-bearing system to base shock excitations using the state-space Newmark scheme and comparisons with experiments［J］.Journal of Sound and Vibration，2006，297:595 －615.

［5］沙云東，聞邦椿，屈 伸.薄壁板在隨機聲載荷作用下的振動響應譜估算［J］.振動與沖擊，2007，26(6):63 －66.

［6］聶旭濤，熊飛嶠.運用統計能量分析法預示空空導彈艙內動力學環境［J］.振動與沖擊，2007，26(4):140 －143.

［7］劉洪英，馮雪梅，馬愛軍.沖擊響應譜控制系統的開發［J］.振動與沖擊，2006，25(6):132 －134.

［8］劉建湖.艦船非接觸水下爆炸動力學的理論與應用［D］.江蘇無錫:中國艦船科學研究中心，2002.

［9］江國和.帶限位器的艦船設備隔離系統沖擊響應研究［D］.上海:上海交通大學，2005.

［10］Lu Y Q，Rui X T.Eigenvalue problem，orthogonal and response of multibody system［A］.Proceedings of the International Conference on Advanced Problems in Vibration Theory and Applications，Beijing:Science Press，2000:711 －714.

［11］Rui X T，Wang G P，Lu Y Q，et al.Transfer matrix method for linear multibody system dynamics［J］.Multibody System Dynamics，2008，19:179 －207.