毛虎平，喬文元，郭保全，王 強(qiáng)，董小瑞

（1.中北大學(xué) 機(jī)械與動力工程學(xué)院，太原 030051； 2.中北大學(xué) 機(jī)電工程學(xué)院，太原 030051）

聚集單元譜元法在承受沖擊載荷結(jié)構(gòu)動態(tài)分析中的應(yīng)用

毛虎平1，喬文元1，郭保全2，王 強(qiáng)1，董小瑞1

（1.中北大學(xué) 機(jī)械與動力工程學(xué)院，太原 030051； 2.中北大學(xué) 機(jī)電工程學(xué)院，太原 030051）

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

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

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

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

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

1 線性動力學(xué)方程的時間譜元法

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

時間譜元法的主要步驟為

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

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

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

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

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

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

1.1 線性結(jié)構(gòu)動態(tài)響應(yīng)方程及其轉(zhuǎn)化形式

線性結(jié)構(gòu)動態(tài)響應(yīng)方程可表達(dá)為

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

1.2 聚集單元劃分

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

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

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

1.3 單元分析

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

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

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

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

通過解式（4）可以獲得GLL點及其權(quán)值。

式中N為插值次數(shù)，權(quán)值可以表示為

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

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

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

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

2 實例分析

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

2.1 標(biāo)準(zhǔn)形式

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

圖2 沖擊載荷

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

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

圖3 等距單元譜元法計算結(jié)果

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

圖4 聚集單元譜元法計算結(jié)果

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

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

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

2.2 線性單自由度吸振器

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

圖7 緩沖器

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

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

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

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

2.3 桿桁架結(jié)構(gòu)

該平面桁架結(jié)構(gòu)包含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負(fù)方向上也作用相同的動態(tài)載荷，b=1 000，ε=0.000 112 5。

圖9 124桿桁架結(jié)構(gòu)

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

圖10 124桿桁架結(jié)構(gòu)關(guān)鍵位置節(jié)點X方向動態(tài)位移響應(yīng)

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

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

2.4 連桿小頭

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

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

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

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

3 結(jié)語

研究結(jié)論如下。

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

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

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

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

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

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

[4]李永強(qiáng)，宦強(qiáng).基于Ansys的沖擊載荷下蜂窩夾芯板的動力學(xué)響應(yīng)[J].東北大學(xué)學(xué)報(自然科學(xué)版)，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]毛虎平，蘇鐵熊，李建軍.基于逐步時間譜元法的結(jié)構(gòu)動態(tài)響應(yīng)仿真[J].中北大學(xué)學(xué)報(自然科學(xué)版)，2013(4)：424-430.

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

[9]張艷崗，毛虎平.動態(tài)應(yīng)力解空間譜元離散的關(guān)鍵時間點識別方法[J].機(jī)械工程學(xué)報，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

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

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