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

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

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

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

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

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

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

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

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

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

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

1 線性動(dòng)力學(xué)方程的時(shí)間譜元法

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

時(shí)間譜元法的主要步驟為

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

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

（3）將每個(gè)時(shí)間單元?jiǎng)澐譃槿舾蓚€(gè)時(shí)間節(jié)點(diǎn)，即正交多項(xiàng)式的零點(diǎn)；

（4）將每個(gè)單元的近似解表示為L(zhǎng)egendre正交多項(xiàng)式的線性組合；

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

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

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

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

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

1.2 聚集單元?jiǎng)澐?/p>

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

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

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

1.3 單元分析

聚集單元?jiǎng)澐趾螅趩卧獏^(qū)間內(nèi)部增加正交多項(xiàng)式的零點(diǎn)，并且可以增加插值點(diǎn)提高近似精度，有限元法中稱(chēng)其為hp法。

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

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

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

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

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

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

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

式中Bu為全局微分矩陣，Bω為全局權(quán)矩陣，為激勵(lì)力的全局形式，Xug為所有狀態(tài)變量在時(shí)間節(jié)點(diǎn)的集合?；?jiǎn)式（9）得

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

2 實(shí)例分析

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

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

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

圖2 沖擊載荷

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

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

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

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

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

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

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

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

2.2 線性單自由度吸振器

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

圖7 緩沖器

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

式（12）自由振動(dòng)的解析解為

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

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

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

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

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

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

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

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

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

2.4 連桿小頭

連桿小頭及其相連的桿身部分，如圖12所示，采用PLANE42單元，右端固定，共有111個(gè)節(jié)點(diǎn)、212個(gè)自由度，在小頭內(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é)點(diǎn)15、42、49、93X方向的位移響應(yīng)。圖13中，將模態(tài)疊加法獲得結(jié)果看作精確解，可以看出在相同單元數(shù)量、尺寸、插值次數(shù)下，聚集單元譜元法更精確，和模態(tài)疊加法保持一致，而等距單元譜元法存在較大誤差。由于位移響應(yīng)本身就很小，因此進(jìn)行了局部放大。

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

3 結(jié)語(yǔ)

研究結(jié)論如下。

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

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

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

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

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

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

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

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

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

國(guó)家自然科學(xué)基金資助項(xiàng)目(51275489)

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