聚集單元譜元法在承受沖擊載荷結構動態(tài)分析中的應用

毛虎平，喬文元，郭保全，王 強，董小瑞

（1.中北大學 機械與動力工程學院，太原 030051； 2.中北大學 機電工程學院，太原 030051）

聚集單元譜元法在承受沖擊載荷結構動態(tài)分析中的應用

毛虎平1，喬文元1，郭保全2，王 強1，董小瑞1

（1.中北大學 機械與動力工程學院，太原 030051； 2.中北大學 機電工程學院，太原 030051）

為了擴展譜元法的適應性，針對承受沖擊載荷的結構動態(tài)響應問題，從譜單元離散方案出發(fā)并根據沖擊載荷的特點，以沖擊載荷最大值點為中心將譜單元尺寸按一定比例等比向兩側擴大，實現單元尺寸與載荷特征相適應，在此基礎上，將動力學方程轉化為1階線性微分方程組，通過Bubnov-Galerkin法獲得離散線性方程組，應用高斯消元法求解。最后，應用此方法對標準形式、線性單自由度吸振器、124桿桁架結構和連桿小頭問題進行計算并與等距譜元法進行比較，說明了此方法的可行性和有效性。

振動與波；聚集單元；譜元法；沖擊載荷；動態(tài)分析

在工程中，機械結構幾乎都要承受動態(tài)載荷，并且多數情況下會承受沖擊載荷，如榔頭敲擊釘子，聲共振無損檢測中敲擊錘敲擊被檢測構件，汽車的碰撞，船舶與橋梁、海洋平臺的碰撞，飛機的墜落，水面艦艇遭受水下非接觸性爆炸的沖擊[1]等，為了評價其動態(tài)特性以及動態(tài)行為有必要進行仿真分析。目前國內外對沖擊載荷作用下結構動態(tài)行為研究較多。比如，起重機械架在起升沖擊載荷作用下的動態(tài)特性分析[2]，沖擊載荷作用下齒輪動態(tài)應力變化研究[3]，沖擊載荷下蜂窩夾芯板的動力響應分析[4]等，很多文獻對各種結構在沖擊載荷的作用下動態(tài)特性及其行為進行了分析研究，然而運用譜元法對沖擊載荷下的結構動態(tài)分析還未見有報道。

譜元法是Patera于1984提出的應用于CFD的一種數值方法，其利用了有限元處理任意結構及其邊界的靈活性和譜方法收斂快的優(yōu)點[5]。譜元法能用較少的單元獲得與其他方法相同的精度，其特點是將每個單元在高斯-勒讓德-洛巴托多項式（Gauss-Legendre-Lobbatto，GLL）的零點處離散，然后進行Lagrange多項式插值。從理論上分析，在正交多項式零點處插值時可獲得最高插值精度[6]。從減小計算規(guī)模角度提出基于逐步時間譜元法的結構動態(tài)響應仿真方法，將仿真時間劃分為很小的時間段，在每一個時間段內劃分單元，每個時間段看做獨立計算部分，前一部分計算結果作為后一部分計算的初始條件，該方案減少了計算時間[7]；在第二類Chebyshev正交多項式零點處，從重心Lagrange插值角度構造了非線性振動問題的離散方案[8]。利用譜單元離散插值精度高的特點，計算了動態(tài)響應優(yōu)化中關鍵時間點[9]。Mohammad H.Kurdi利用時間譜元法求解了簡單的質量彈簧阻尼系統，在此基礎上對單自由度吸振器和單自由度微型控制器進行了優(yōu)化[10]。

本文從結構動力學控制方程出發(fā)，利用有限元進行空間離散，用譜元法進行時間離散，其中針對沖擊載荷時間短、變化大的特點，將譜單元的離散進行聚集處理，彌補等距單元誤差大的缺陷。

1 線性動力學方程的時間譜元法

譜元法分為時間譜元法、空間譜元法和時空混合譜元法。對于動力學方程而言，可以采用三種方法中的任何一種，文中利用譜元法離散時間，采用有限元法離散空間。

時間譜元法的主要步驟為

（1）將動力學方程轉化為一階線性微分方程組，然后通過布巴諾夫-伽遼金法（Bubnov-Galerkin）等效為積分形式，代入單元積分表達式，獲得時間單元譜元方程；

（2）根據沖擊載荷的特點，將仿真時間劃分為聚集型單元，即在載荷突變的位置單元尺寸比較小，在載荷平坦的位置單元尺寸比較大，這可通過單元最大尺寸和最小尺寸的比值來控制；

（3）將每個時間單元劃分為若干個時間節(jié)點，即正交多項式的零點；

（4）將每個單元的近似解表示為Legendre正交多項式的線性組合；

（5）利用連接矩陣將所有單元的近似解組裝成總體譜元方程；

（6）求解總體譜元方程，獲得全局位移近似解和速度近似解，可以通過動力學微分方程得到加速度的解，對于結構來說，可以通過應力和位移的關系獲得節(jié)點應力解。

1.1 線性結構動態(tài)響應方程及其轉化形式

線性結構動態(tài)響應方程可表達為

其中M為質量矩陣，C為阻尼矩陣，K為剛度矩陣，F為動態(tài)載荷向量，x為位移向量，為加速度向量，為速度向量。初始條件為x(0)=b0，其中M、C、K不隨時間變化；是時間的函數，F是任意時間函數

1.2 聚集單元劃分

如圖1所示，將仿真時間t∈[t0,tn]分割為n個互不相交的單元，即每個單元配置若干個點。在劃分單元時，以沖擊載荷最大值點為中心，中心處單元尺寸最小，越往兩邊尺寸越大，如圖1所示。

圖1 聚集單元譜元法離散時間及對應的沖擊載荷

每個單元配置Chebyshev或Legendre正交多項式的零點，其點數可以相同，也可不同，文中采用配置相同的零點。這樣可以避免由于沖擊載荷的局部突變性帶來的求解誤差。具體實施時，通過設來實現，圖1中ratio=0.1。圖1 (a)為聚集單元示意圖，圖1(b)表示30個聚集單元對應的沖擊載荷，其中沖擊載荷可以有3個參數控制，即沖擊載荷最大值、寬度和載荷中心。

1.3 單元分析

聚集單元劃分后，在單元區(qū)間內部增加正交多項式的零點，并且可以增加插值點提高近似精度，有限元法中稱其為hp法。

在聚集單元譜元法中，通過在單元內特殊位置配置插值GLL點，其基函數表示為正交多項式的組合，構成單元內各點的形函數，這樣可以在有限個點上插值達到譜方法的收斂速度。

在時間譜元法中一般采用兩種正交多項式，即Chebyshev和Legendre正交多項式，文中采用Legendre正交多項式[6]。

由于單元端點不是Legendre正交多項式的零點，因此加入Lobbato多項式保證了單元端點包含在零點中。Lobbato多項式滿足正交特性

通過解式（4）可以獲得GLL點及其權值。

式中N為插值次數，權值可以表示為

將狀態(tài)變量在仿真時間內按圖1所示方案并結合沖擊載荷離散，表示為m次Lagrange多項式

集成是把所有單元的譜元方程按照離散中的順序組合起來，獲得總體譜元方程。有Nv狀態(tài)方程的系統，通過耦合矩陣的方陣)的張量叉乘得到全部狀態(tài)變量的全局組裝

式中Bu為全局微分矩陣，Bω為全局權矩陣，為激勵力的全局形式，Xug為所有狀態(tài)變量在時間節(jié)點的集合。化簡式（9）得

式中G為時間段的全局線性矩陣，式（10）是線性方程組，通過高斯消元法求解。

2 實例分析

為了驗證聚集單元譜元法的優(yōu)勢，對4個不同實例進行分析，并與等距單元譜元法比較。

2.1 標準形式

對于突變載荷，常規(guī)的等長度單元譜元法不能獲得好的收斂效果。沖擊載荷如圖2所示。

圖2 沖擊載荷

無論多么復雜的系統，動態(tài)響應方程最終可以轉化為形如式（11）所示的標準形式。其中，載荷在0.5 s時刻，產生很大的沖擊，沖擊的幅值可以由系數決定，沖擊力的作用時間可以通過調整ε來改變，本例中ε=0.000 112 5。解析解可以通過積分因子法并借用Matlab中的erf函數獲得

圖3由等距單元譜元法獲得，很明顯誤差很大，特別在載荷突變處，Ne=30時，誤差極大，而當Ne=100時，誤差減小了很多，但還是較大。

圖3 等距單元譜元法計算結果

圖4由聚集單元譜元法獲得，在單元數和插值次數相等的前提下已經消除了誤差。

圖4 聚集單元譜元法計算結果

圖5和圖6更能說明聚集單元譜元法的優(yōu)勢，圖5將等距單元譜元法與聚集單元譜元法比較，發(fā)現在50個單元處，單元插值次數為18時前者為（50，0.02506）；而單元插值次數為8時，后者為（50，1.382×10-10），有非常大的優(yōu)勢；圖6將聚集單元譜元法與等距單元譜元法比較，發(fā)現在80個單元處，單元插值次數分別為5、10、18時前者為（80，3×10-9）、（80，3.628× 10-12）、（80，2.668×10-12），而后者在80個單元處，單元插值次數為18次時為（80，0.011 24），可以看出聚集單元譜元法求解沖擊載荷響應時有很大優(yōu)勢。

圖5 聚集單元與等距單元比較

圖6 等距單元與聚集單元比較

2.2 線性單自由度吸振器

圖7為包含固定質量m=1 kg、彈簧剛度系數和阻尼器系數分別為k=0.9、c=0.9的線性單自由度系統。

圖7 緩沖器

當t=0時，系統以速度v=1 m/s撞擊在一個固定阻礙物上，同時受沖擊載荷作用，ε=0.011 25。運動方程為

式（12）自由振動的解析解為

線性單自由度吸振器雖然既有初始速度，又同時受沖擊載荷，但沖擊載荷的中心在6 s處，寬度由ε=0.011 25來控制而且很窄，因此在4 s附近還是屬于自由振動，可以與自由振動解析解比較，如圖8所示，從圖中可以看出單元數為12，單元插值次數為12時，聚集單元譜元法明顯比等距單元譜元法精度高，說明本文方法具有優(yōu)勢。

圖8 線性單自由度吸振器動態(tài)位移響應

2.3 桿桁架結構

該平面桁架結構包含49個鉸鏈，94個自由度，見圖9。彈性模量E=207 GPa，泊松比v=0.3，密度ρ=7 850 kg/m3，桿的截面積0.645×10-4m2。在節(jié)點1、20、19、18、17、16、15的x正方向上作用相同的動態(tài)載荷，在節(jié)點1、2、3、4、5的y負方向上也作用相同的動態(tài)載荷，b=1 000，ε=0.000 112 5。

圖9 124桿桁架結構

從工程角度，將節(jié)點1、2、3、4、5、15、16、17、18、19、20節(jié)點作為重點考察的位置，稱之為關鍵位置，同時也是沖擊載荷作用的位置。在沖擊載荷作用下，其位移響應如圖10所示，將其中節(jié)點1的X方向位移單獨考慮，并局部放大如圖11所示。

圖10 124桿桁架結構關鍵位置節(jié)點X方向動態(tài)位移響應

圖11 124桿桁架結構節(jié)點1X方向動態(tài)位移響應

可以明顯看出等距單元譜元法計算結果有部分點偏離精確值，而聚集單元譜元法沒有這種現象，實際上關鍵位置的每一個點都出現類似現象，限于篇幅，沒有全部列出。

2.4 連桿小頭

連桿小頭及其相連的桿身部分，如圖12所示，采用PLANE42單元，右端固定，共有111個節(jié)點、212個自由度，在小頭內圓內中間受水平向右的沖擊載荷，載荷如圖2所示，b=106，ε=0.000 112 5，連桿材料的彈性模量E= 207 GPa，泊松比v=0.3，密度ρ=7 850 kg/m3。

圖12 連桿小頭與桿身相連部分

連桿小頭圓角處為關鍵位置，考察節(jié)點15、42、49、93X方向的位移響應。圖13中，將模態(tài)疊加法獲得結果看作精確解，可以看出在相同單元數量、尺寸、插值次數下，聚集單元譜元法更精確，和模態(tài)疊加法保持一致，而等距單元譜元法存在較大誤差。由于位移響應本身就很小，因此進行了局部放大。

圖13 連桿4個關鍵位置節(jié)點X方向動態(tài)位移響應

3 結語

研究結論如下。

（1）針對沖擊載荷的特點，并結合譜元法精度高的優(yōu)點，采用設定載荷中心單元尺寸與兩側單元尺寸的比值來適應載荷的突變性，使單元的離散與沖擊載荷變化相適應；

（2）在單元數量、尺寸、插值次數相同的前提下，求解沖擊載荷動態(tài)響應，聚集單元譜元法的精度比等距單元譜元法高；

（3）從動態(tài)問題的標準形式，到線性單自由度吸振器、124桿平面桁架，再到連桿小頭問題的沖擊載荷動態(tài)響應分析，都說明了本文方法的可行性和有效性。這為進一步研究沖擊載荷作用下的結構動態(tài)分析提供參考。

[1]肖鋒，諶勇，馬超，等.橡膠蜂窩覆蓋層水下爆炸響應及抗沖擊性能[J].噪聲與振動控制，2013，33(4)：44-49.

[2]王貢獻，沈榮瀛.起重機臂架在起升沖擊載荷作用下動態(tài)特性研究[J].機械強度，2005，27(5)：561-566.

[3]唐進元，彭方進，黃云飛.沖擊載荷下的齒輪動應力變化規(guī)律數值分析[J].振動與沖擊，2009，28(08)：138-143.

[4]李永強，宦強.基于Ansys的沖擊載荷下蜂窩夾芯板的動力學響應[J].東北大學學報(自然科學版)，2015：36(6)：858-862.

[5]PATERA A T.A spectral element method for fluid dynamics:Laminarflow in achannelexpansion[J].Journal of Computational Physics,1984,54(3):468-488.

[6]POZRIKIDIS C.Introduction to finite and spectral element methods using Matlab[M].CRC Press,2005.

[7]毛虎平，蘇鐵熊，李建軍.基于逐步時間譜元法的結構動態(tài)響應仿真[J].中北大學學報(自然科學版)，2013(4)：424-430.

[8]毛虎平，王偉能，續(xù)彥芳，等.非線性振動分析的切比雪夫譜元法[J].噪聲與振動控制，2015，35(1)：73-77.

[9]張艷崗，毛虎平.動態(tài)應力解空間譜元離散的關鍵時間點識別方法[J].機械工程學報，2014，50(5)：82-87.

[10]KURDI M H,BERAN P S.Optimization of dynamic response using a monolithic-time formulation[J].Structural and Multidisciplinary Optimization,2009, 39(1):83-104.

Application of SEM based onAccumulation Elements to DynamicAnalysis of the Structures Subjected to Impact Loads

MAO Hu-ping1,QIAO Wen-yuan1,GUO Bao-quan2, WANG Qiang1,DONG Xiao-rui1
(1.College of Mechanical and Power Engineering,North University of China,Taiyuan 030051,China; 2.College of Mechatronic Engineering,North University of China,Taiyuan 030051,China)

The spectral element method(SEM)is used to analyze the dynamic responses of the structures subjected to impact loading.Starting with the SEM discrete scheme and according to the feature of the impact loads,each spectral element is enlarged proportionally from the middle point of the impact load to both sides of the element to realize the suitability between the element size and the loading feature.On this basis,the dynamic equation is transformed to a set of the first-order linear differential equations,and the corresponding discrete linear algebraic equations are obtained by Bubnov-Galerkin method and solved by Gauss elimination method.Finally,four examples of a typical absorber,a linear single DOF absorber,a 124-member plane truss and the small end of a link are calculated with this method and the results are compared with those of the equal-spaced SEM.The feasibility and validity of this method are verified.

vibration and wave;accumulation element;spectral element method(SEM);impact load;dynamic analysis

TH122

：A

：10.3969/j.issn.1006-1335.2016.06.009

1006-1355(2016)06-0045-06

2016-07-04

國家自然科學基金資助項目(51275489)

毛虎平（1974-），男，博士，副教授，碩士生導師，主要從事振動理論與工程數值分析方法、結構動態(tài)響應優(yōu)化方法研究。E-mail:maohp@nuc.edu.cn