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

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

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

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

譜元法分為時間譜元法、空間譜元法和時間-空間耦合譜元法。對于動力學(xué)方程而言，可以采用其中任何一種方法。下面以時間譜元法為例進行介紹。

時間譜元法的主要步驟：①將動力學(xué)方程轉(zhuǎn)化為一階線性微分方程組，然后通過Bubnov-Galerkin法等效為積分形式，代入單元積分表達式，獲得時間單元譜元方程；②根據(jù)沖擊載荷的特點，將仿真時間劃分為聚集型單元，即在載荷突變的位置單元尺寸比較小，在載荷平坦的位置單元尺寸比較大，它可通過單元最大尺寸和最小尺寸的比值來控制；③將每個時間單元劃分為若干個時間節(jié)點，即正交多項式的零點；④將每個單元的近似解表示為Legendre正交多項式的線性組合；⑤利用連接矩陣將所有單元的近似解集成總體譜元方程；⑥求解總體譜元方程，獲得全局位移近似解和速度近似解，可以通過動力學(xué)微分方程得到加速度的解，對于結(jié)構(gòu)來說，可以通過應(yīng)力和位移的關(guān)系獲得節(jié)點應(yīng)力解。

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

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

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

為了采用時間譜元法，設(shè)x1=，x2=x，式（2.38）可轉(zhuǎn)化為一階線性常微分方程組［見式（2.39）］，此方程組與式（2.38）同解。

2.3.2 聚集單元劃分

如圖2.5所示，將仿真時間t∈［t 0，tn］分割為n個互不相交的單元，即［t0，t1］，［t1，t2］，［t 2，t3］，…，［t n-1，tn］，每個單元配置若干個點。在劃分單元時，以沖擊載荷最大值點為中心，中心處單元尺寸最小，越往兩邊單元尺寸越大。每個單元配置Chebyshev或Legendre正交多項式的零點，其點數(shù)可以相同，也可以不同，本章采用相同的零點配置。這樣可以避免沖擊載荷的局部突變帶來的求解誤差。在具體實施時，可以設(shè)（圖2.5中ratio=0.1）。圖2.5（a）所示為聚集單元示意，圖2.5（b）所示為30個聚集單元對應(yīng)的沖擊載荷，其中，沖擊載荷由3個參數(shù)控制，即沖擊載荷最大值、寬度和載荷中心。

2.3.3 單元分析

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

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

在時間譜元法中一般采用兩種正交多項式，即Chebyshev和Legendre正交多項式，本章采用Legendre正交多項式［86］。

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

式中，δij是Kronecker δ函數(shù)；

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

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

式中，ξk為LoN-1（ξ）的零點，k=1，2，…，N-1。

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

式中，為j單元的k次Lagrange多項式；ξk為定義在［-1，1］上的GLL點；是單元j上未知節(jié)點在GLL點的值。然后通過Bubnov-Galerkin法可得

用矩陣形式表示每個單元的離散：

式中，，Iω為ξ的一般函數(shù)，

2.3.4 集成總體譜元方程及求解

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

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

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