基于減基法的有阻尼結構諧響應快速分析方法

黃 芬 韓 旭 黃永輝

1.湖南大學汽車車身先進設計制造國家重點實驗室,長沙,410082 2.中聯重工科技股份有限公司,長沙,410013

0 引言

結構諧響應分析是結構分析很重要的一方面,它可以確保一個給定的結構能經受住各種不同頻率的正弦載荷,有助于結構設計人員掌握結構的持續動力特性,并在必要時避免或利用共振。結構諧響應分析需要求解大型多自由度系統方程,這個循環求解的過程非常耗時。目前有許多降階方法可用來解決這類問題,如Guyan靜力凝聚法[1-2]、動力子結構法[3]等,但這些方法在提高效率和保持原模型的物理特性上仍存在一定的局限性,即使是降階后,求解仍比較耗時。

減基法[4-5]是20世紀70年代提出的一種降階方法。其基本思想是,當系統由多個參數描述時,這些參數的不同組合會使系統方程有不同的解,而系統在新參數下的解可以由事先設計的樣本參數組所對應解的線性組合來得到。近年來,不少學者對該方法進行了研究。國外,Maday等[6]從理論上證明了減基法的一致指數收斂性質;Rovas[7]把減基法應用于不同種類偏微分方程的求解中;國內,減基法在結構動力學中的應用研究處于起步階段,文獻[8]引入減基法,在波數域內快速分析了復合材料層合板的瞬態響應問題。李永紅等[9]將減基法用于結構特征值分析與參數優化中,求解了無阻尼結構的固有頻率與振型。黃永輝等[10]通過減基法在實數范圍內分析了無阻尼結構的諧響應問題。然而,在目前可查資料中,尚未見減基法用于有阻尼結構諧響應分析的報導。

本文在有阻尼的線彈性結構諧響應分析中引入減基法,在離散化的頻域中進行對數預采樣[8],求出樣本點的解向量,得到一系列的復數解向量。而無阻尼的結構諧響應分析得到的樣本點的解向量是一系列的實數解向量,后續工作都是在實數域中進行的。有阻尼的情況要考慮直接用復數解向量或是向量的模來構造減基空間。復數向量包含幅值和相位兩個屬性,用模來構造減基空間,把復數域的工作轉換到實數域進行,會丟失向量的相位屬性,因而考慮直接利用復數解向量,通過貪婪算法[6]構造減基空間,把系統的剛度矩陣、質量矩陣、阻尼矩陣以及載荷列向量分別投影到減基空間進行降階,得到減縮系統。在減基空間求解原問題的減縮解,并把減縮解還原到原空間得到原問題的近似解。通過誤差和計算耗時的對比,驗證了這種方法在有阻尼的結構諧響應分析中的可靠性和有效性。

1 有阻尼結構諧響應分析的有限元形式

諧響應分析是用來計算結構在簡諧載荷激勵下響應的方法。在諧響應分析中,激勵外載是已知的,可以是力也可以是位移、速度或加速度。

在諧響應分析中,激勵載荷在本質上為正弦載荷,在最簡單的情況下,這種載荷可通過特定頻率下的幅值來定義。簡諧響應與載荷以相同的頻率出現,由于系統阻尼的影響,響應在時間上可能移位,響應移位也稱相位移位,因為載荷峰值與響應峰值不是同時出現的。

在諧響應分析中,一般有兩種不同的數值方法可供選擇：直接法和模態法。直接法根據外載荷頻率求解耦合的運動方程;模態法利用結構的振型縮減和解耦運動方程,對各個模態響應進行疊加得到一特定外載荷頻率的解。

簡諧載荷作用下的有阻尼結構的運動方程為

其中,F0(t)為簡諧激勵載荷,F0(t)=F(ω)eiωt;F為載荷幅值;ω為激勵載荷的頻率,當載荷幅值不變時,F與 ω無關,可寫成F0(t)=F eiωt的形式。K、M、C分別為結構的整體剛度矩陣、質量矩陣和分別為結構的整體位移向量、速度向量、加速度向量。對于簡諧運動,假設一個簡諧形式的解：

式中,d0(ω)為復位移向量。

將式(2)代入式(1)并化簡得

當考慮阻尼時,式(3)代表復系數方程系統。當激勵頻率范圍較大時,式(3)的求解次數也隨著增多,反復求解工作量大。另一方面,當系統自由度n很大時,求式(3)非常耗時,需要考慮計算成本問題。

在實際問題中,要精確地計算結構受到的阻尼是很困難的。通常從整體上考慮阻尼的影響,近似地估計阻尼力做功消耗的能量。一般有限元程序中,假定整體阻尼矩陣C是整體剛度矩陣K與整體質量矩陣M的線性組合,成為Rayleigh阻尼(或比例阻尼),其表達式為

其中,Ω為結構固有頻率;ξ為阻尼比;比例系數α、β可通過測試結構自由振動的衰減率經過換算得到。當阻尼比ξ確定時,α、β也可通過式(5)得到。

2 有阻尼結構的減基法諧響應分析

減基法用于結構諧響應分析的目的是縮減剛度矩陣、質量矩陣和阻尼矩陣的階數,同時保證減縮后的系統能很好地近似原系統,以提高求解效率。計算過程分為離線和在線兩階段。

2.1 離線階段

結構諧響應分析的主要變化參數是載荷頻率ω。本文將頻域離散成m個頻率點,對頻率點進行預采樣,目前預采樣方法有對數采樣法、等間隔均勻采樣法、Chebychev采樣方法等。本文用對數采樣方法采取N個樣本點,采樣公式如下：

式中,λ為常數;μj為采樣參數的值;μmax為采樣參數的最大值。

N＜m時,初步構造頻率的樣本空間如下：

在每一個樣本點求式(3),由于本文考慮的是有阻尼的系統,式(3)中存在復數項,所以得到的N個解向量也是復向量,樣本點的復數解向量構成解空間如下：

將d0(ω)分離成與變化參數有關部分 α(ω)和無關部分D,則有

那么式(3)可寫成

式中,α為原系統在減基空間的減縮解。

兩邊乘以DT得

當D為n ×N1階正交矩陣時,稱K N1、C N1、M N1、F N1分別為K、C、M、F關于D的正交投影,且K、M、C、F的階數由n×n變成N 1×N1,當K N1=DTKD,CN1=DTCD,MN1=DTMD,FN1=DTF,N1＜n時,矩陣得到降階。

為了得到正交矩陣D,構造減基空間,本文引入貪婪算法,對解空間WN進行自適應訓練。這個過程中,與無阻尼情況的區別是,樣本點的解向量都是復數向量,構成的空間也是復空間,貪婪算法是在復數域進行的。為了保證解向量的正交性,同時避免減縮矩陣的病態問題,每將一個解向量加入減基空間,就要采用施密特正交化方法對減基空間的向量進行正交歸一化。貪婪算法的具體步驟如下：

(1)從預采樣點的解空間中選擇一個解向量加入減基空間。

(2)當構造減基空間的向量個數大于1時,用施密特正交化方法把解向量正交歸一化,構造新的減基空間,并得到正交矩陣D。

(3)通過伽遼金映射將與參數 ω無關的K、M、C、F矩陣投影到減基空間,得到縮減后的矩陣和系統控制方程,在預采樣點進行求解,得到 N個解向量。

4.假設該化妝品使用的稅率如下：進口關稅稅率=t，消費稅率=c,增值稅率=ν。根據《關于跨境電子商務零售進口稅收政策的通知》，計算可得：

(4)利用2范數,求N個樣本點的減基誤差和投影誤差。減基誤差向量e r的第j個分量定義為

投影誤差向量e p的第j個分量定義為

式中,d0j為直接求解的第j個樣本點的解向量;αj為減基空間中的第j個樣本點的解向量。

(5)判斷最大減基誤差是否小于給定的誤差限ε,若滿足,停止貪婪算法;若不滿足,將最大投影誤差對應的預采樣點的解向量放入減基空間,重復步驟(2)～步驟(4)。

通過貪婪算法構造的減基空間能很好地近似原來的空間,可表示為

式中,ηi(i=1,2,…,N1)為構成減基空間的樣本點的解向量。

且一般情況下N1＜N。

2.2 在線階段

利用正交矩陣D將K、M、C、F降階得減縮矩陣 K N1、C N1、M N1、F N1。

由于減基空間中的各個向量線性無關,故這些向量可以作為任意參數下解空間中N 1維子空間的基,則原系統的解向量可以用減基空間解向量的線性疊加來近似：

3 減基法與直接法之間的誤差

為了驗證減基法的有效性,需要定義減基法求解與直接求解之間的誤差,利用歐幾里德范數定義相對誤差：

其中,d0k為直接求解式(3)得到的解向量,根據數值分析中方程組近似解可靠性判別法[11],可利用系數矩陣條件數以及殘余向量來判斷減基法得到的近似解是否可靠。令

當

時,近似解是可靠的。其中,Cond(*)為求矩陣A條件數的函數。

4 算例

算例一 以一L形振動試驗夾具為例,材料常數如下：彈性模量 E=70GPa,泊松比μ=0.33,密度ρ=2.8×103kg/m3,夾具板厚度0.01m。該夾具通過4個螺栓連接在振動平臺臺面上,其幾何結構尺寸及螺栓位置如圖1所示,有限元模型如圖2所示,共有791個節點、2373個自由度,在747節點y方向施加簡諧載荷,幅值1000N,頻率3500～4500Hz,系統受瑞利阻尼作用,阻尼比ξ=0.05。

將頻域離散成200個子步,用對數預采樣采集31個頻率點,給定的減基誤差限ε1=1×10-5,經過貪婪算法,最終選出7個點即可滿足精度要求。減基法把原來2373×2373的系數矩陣縮減為7×7的矩陣,大幅度提高了求解效率。計算用MATLAB程序實現,在CPU主頻為2.01GHz、內存為1GB的電腦上運行,直接求解與減基法方法求解時間對比見表1。減基法求解所耗時間幾乎是直接求解的1/5,效率明顯提高。

圖 1 夾具幾何尺寸(單位：mm)

圖2 夾具有限元模型

貪婪算法過程中,減基誤差、投影誤差分別如圖3、圖4所示,減基誤差從第1代到第2代迅速減小,隨后逐漸收斂于誤差限10-5。投影誤差變化趨勢與減基誤差基本相同,第3代后緩慢收斂,繼續增加基數量對誤差影響不大,以減基誤差判斷,當基數量達到7時,誤差小于給定誤差限。

圖3 減基誤差隨著基向量個數增加的變化情況

相對誤差如圖5所示,其中最大相對誤差為3.73×10-11,滿足式(18),圖6給出了減基法求解與直接求解得到的響應幅值,找到一頻率在4000Hz附近的共振點,兩種所得結果一致,說明減基法能得到高精度的解,且所耗時間比直接求解所耗時間少。

圖4 投影誤差隨著基向量個數增加的變化情況

圖5 減基法求解與直接求解之間的相對誤差

圖6 減基法與原方法的響應幅值比較(701節點y方向)

算例二 考慮如圖7所示的板梁支架結構,梁底端固支,彈性模量E=200GPa,泊松比μ=0.3,密度 ρ=7.8×103kg/m3,板長 0.8m,寬0.4m,厚 0.01m。梁截面為實心,寬0.03m,高0.04m。支架的每段支撐高0.4m。有限元模型如圖7所示,共有405個節點、2430個自由度,分別在支架的219節點、270節點、347節點加載z方向同頻率簡諧載荷,幅值分別為1000N、750N、500N,頻率為0～300Hz,同樣受瑞利阻尼作用,阻尼比ξ=0.06。

與算例一同理,在離散為300個子域的頻域中采取16個點,貪婪算法選出10個點即可滿足精度要求,本例減基誤差限ε=1×10-3,兩種方法耗時如表2所示。減基誤差、投影誤差的變化如圖8、圖9所示,趨勢基本相同,從第 1代到第7代迅速減小,隨后逐漸收斂于誤差限10-2。當基數量達到11時,誤差小于給定誤差限,此時,構造的減基空間滿足精度要求,減縮系統的剛度矩陣、質量矩陣、阻尼矩陣階次由 2430×2430降到11×11。

圖7 支架結構有限元模型

表2 減基法求解和直接求解的時間對比

圖8 減基誤差隨著基向量個數增加的變化情況

圖9 投影誤差隨著基向量個數增加的變化情況

最后的相對誤差如圖10所示,其中最大誤差為1.01×10-5,滿足式(18),同時由表2知,減基法耗時只是直接求解耗時的1/7。圖11給出了減基法與直接法求解的響應幅值,實線為有限元法直接求解的結果,點線為減基法求解結果,由圖11可見減基法與原有限元法的結果一致,通過曲線找到共振點頻率約為80Hz。

以上兩個算例說明在有阻尼的線彈性結構諧響應分析中,利用減基法求解可以明顯提高效率,并保證解的精度,最大誤差均滿足誤差限,說明減基法和原方法的結果是可靠的,也說明了該方法的有效性。

5 結束語

圖10 減基法求解與直接求解之間的相對誤差

圖11 減基法與原方法的響應幅值比較(86節點 x方向)

本文在有限元的基礎上,把減基法用于有阻尼的線彈性結構諧響應分析。分析過程中,把簡諧載荷頻率離散,通過對數預采樣和貪婪算法構造減基空間,把原系統映射到減基空間進行降階,在很大程度上縮減了系統剛度矩陣、質量矩陣、阻尼矩陣的階數,從而節省了求解資源,有效提高效率。本文算例表明,減基法在解決無阻尼的結構諧響應分析的基礎上,同樣適用于有阻尼的結構諧響應分析。

[1] Gu J.Efficient Model Reduction Methods for Structural Dynamics Analyses[D].Michigan：The University of Michigan,2000.

[2] Guyan R J.Reduction of Stiffness and Mass Matrices[J].AIAA Joumal,1965,3(2)：380-381.

[3] Wilson E L,Yuan M W,Dickens J M.Dynamic Analysis by Direct Superposition of Ritz Vectors[J].Earthquake Engineering Structural Dynamics,2006,10(6)：813-821.

[4] Noor A K,Peters JM.Reduced Basis Technique for Nonlinear Analysis of Structures[J].AIAA Journal,1980,18(4)：455-462.

[5] Veroy K,Patera A T.Certified Real-time Solution of the Parameterized Steady Incompressible Navier-Stokes Equations：Rigorous Reduced-basis a Posteriori Error Bounds[J].International Journal for Numerical Method in Fluids,2005,47(8/9)：773-788.

[6] Maday Y,Patera A T,Turinici G.A Priori Convergence Theory for Reduced–basis Approximations of Single-parameter Elliptic Partial Differential E-quations[J].Journal of Scientific Computing,2002,17(1/4)：437-446.

[7] Rovas DV.Reduced-basis Output Bound Methods for Parametrized Partial Differential Equations[D].Cambridge：Massachusetts Institute of Technology,2002.

[8] 黃永輝,韓旭,冉承新.基于減基法的層合板瞬態響應快速分析方法[J].力學學報,2008,40(2)：255-260.

[9] 李永紅,李光耀.Reduced-basis Method在特征值問題中的應用[J].機械制造,2008,46(525)：5-9.

[10] 黃永輝,韓旭,黃芬.基于減基法的結構諧響應快速分析方法[J].振動與沖擊,2009,28(7)：61-64.

[11] Gerald C F,Wheatley PO.應用數值分析[M].北京：機械工業出版社,2006.

[12] 馬超,呂震宙.結構可靠性分析的支持向量機分類迭代算法[J].中國機械工程,2007,18(7)：816-819.

[13] 許兆堂,朱如鵬.阻尼分段階梯傳動軸主共振的分析[J].中國機械工程,2006,17(7)：685-690.