基于減基法的有阻尼結(jié)構(gòu)諧響應(yīng)快速分析方法

黃 芬 韓 旭 黃永輝

1.湖南大學(xué)汽車車身先進(jìn)設(shè)計(jì)制造國家重點(diǎn)實(shí)驗(yàn)室,長(zhǎng)沙,410082 2.中聯(lián)重工科技股份有限公司,長(zhǎng)沙,410013

0 引言

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

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

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

1 有阻尼結(jié)構(gòu)諧響應(yīng)分析的有限元形式

諧響應(yīng)分析是用來計(jì)算結(jié)構(gòu)在簡(jiǎn)諧載荷激勵(lì)下響應(yīng)的方法。在諧響應(yīng)分析中,激勵(lì)外載是已知的,可以是力也可以是位移、速度或加速度。

在諧響應(yīng)分析中,激勵(lì)載荷在本質(zhì)上為正弦載荷,在最簡(jiǎn)單的情況下,這種載荷可通過特定頻率下的幅值來定義。簡(jiǎn)諧響應(yīng)與載荷以相同的頻率出現(xiàn),由于系統(tǒng)阻尼的影響,響應(yīng)在時(shí)間上可能移位,響應(yīng)移位也稱相位移位,因?yàn)檩d荷峰值與響應(yīng)峰值不是同時(shí)出現(xiàn)的。

在諧響應(yīng)分析中,一般有兩種不同的數(shù)值方法可供選擇：直接法和模態(tài)法。直接法根據(jù)外載荷頻率求解耦合的運(yùn)動(dòng)方程;模態(tài)法利用結(jié)構(gòu)的振型縮減和解耦運(yùn)動(dòng)方程,對(duì)各個(gè)模態(tài)響應(yīng)進(jìn)行疊加得到一特定外載荷頻率的解。

簡(jiǎn)諧載荷作用下的有阻尼結(jié)構(gòu)的運(yùn)動(dòng)方程為

其中,F0(t)為簡(jiǎn)諧激勵(lì)載荷,F0(t)=F(ω)eiωt;F為載荷幅值;ω為激勵(lì)載荷的頻率,當(dāng)載荷幅值不變時(shí),F與 ω?zé)o關(guān),可寫成F0(t)=F eiωt的形式。K、M、C分別為結(jié)構(gòu)的整體剛度矩陣、質(zhì)量矩陣和分別為結(jié)構(gòu)的整體位移向量、速度向量、加速度向量。對(duì)于簡(jiǎn)諧運(yùn)動(dòng),假設(shè)一個(gè)簡(jiǎn)諧形式的解：

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

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

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

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

其中,Ω為結(jié)構(gòu)固有頻率;ξ為阻尼比;比例系數(shù)α、β可通過測(cè)試結(jié)構(gòu)自由振動(dòng)的衰減率經(jīng)過換算得到。當(dāng)阻尼比ξ確定時(shí),α、β也可通過式(5)得到。

2 有阻尼結(jié)構(gòu)的減基法諧響應(yīng)分析

減基法用于結(jié)構(gòu)諧響應(yīng)分析的目的是縮減剛度矩陣、質(zhì)量矩陣和阻尼矩陣的階數(shù),同時(shí)保證減縮后的系統(tǒng)能很好地近似原系統(tǒng),以提高求解效率。計(jì)算過程分為離線和在線兩階段。

2.1 離線階段

結(jié)構(gòu)諧響應(yīng)分析的主要變化參數(shù)是載荷頻率ω。本文將頻域離散成m個(gè)頻率點(diǎn),對(duì)頻率點(diǎn)進(jìn)行預(yù)采樣,目前預(yù)采樣方法有對(duì)數(shù)采樣法、等間隔均勻采樣法、Chebychev采樣方法等。本文用對(duì)數(shù)采樣方法采取N個(gè)樣本點(diǎn),采樣公式如下：

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

N＜m時(shí),初步構(gòu)造頻率的樣本空間如下：

在每一個(gè)樣本點(diǎn)求式(3),由于本文考慮的是有阻尼的系統(tǒng),式(3)中存在復(fù)數(shù)項(xiàng),所以得到的N個(gè)解向量也是復(fù)向量,樣本點(diǎn)的復(fù)數(shù)解向量構(gòu)成解空間如下：

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

那么式(3)可寫成

式中,α為原系統(tǒng)在減基空間的減縮解。

兩邊乘以DT得

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

為了得到正交矩陣D,構(gòu)造減基空間,本文引入貪婪算法,對(duì)解空間WN進(jìn)行自適應(yīng)訓(xùn)練。這個(gè)過程中,與無阻尼情況的區(qū)別是,樣本點(diǎn)的解向量都是復(fù)數(shù)向量,構(gòu)成的空間也是復(fù)空間,貪婪算法是在復(fù)數(shù)域進(jìn)行的。為了保證解向量的正交性,同時(shí)避免減縮矩陣的病態(tài)問題,每將一個(gè)解向量加入減基空間,就要采用施密特正交化方法對(duì)減基空間的向量進(jìn)行正交歸一化。貪婪算法的具體步驟如下：

(1)從預(yù)采樣點(diǎn)的解空間中選擇一個(gè)解向量加入減基空間。

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

(3)通過伽遼金映射將與參數(shù) ω?zé)o關(guān)的K、M、C、F矩陣投影到減基空間,得到縮減后的矩陣和系統(tǒng)控制方程,在預(yù)采樣點(diǎn)進(jìn)行求解,得到 N個(gè)解向量。

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

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

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

式中,d0j為直接求解的第j個(gè)樣本點(diǎn)的解向量;αj為減基空間中的第j個(gè)樣本點(diǎn)的解向量。

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

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

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

且一般情況下N1＜N。

2.2 在線階段

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

由于減基空間中的各個(gè)向量線性無關(guān),故這些向量可以作為任意參數(shù)下解空間中N 1維子空間的基,則原系統(tǒng)的解向量可以用減基空間解向量的線性疊加來近似：

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

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

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

當(dāng)

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

4 算例

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

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

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

圖2 夾具有限元模型

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

圖3 減基誤差隨著基向量個(gè)數(shù)增加的變化情況

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

圖4 投影誤差隨著基向量個(gè)數(shù)增加的變化情況

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

圖6 減基法與原方法的響應(yīng)幅值比較(701節(jié)點(diǎn)y方向)

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

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

圖7 支架結(jié)構(gòu)有限元模型

表2 減基法求解和直接求解的時(shí)間對(duì)比

圖8 減基誤差隨著基向量個(gè)數(shù)增加的變化情況

圖9 投影誤差隨著基向量個(gè)數(shù)增加的變化情況

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

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

5 結(jié)束語

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

圖11 減基法與原方法的響應(yīng)幅值比較(86節(jié)點(diǎn) x方向)

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

[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] 黃永輝,韓旭,冉承新.基于減基法的層合板瞬態(tài)響應(yīng)快速分析方法[J].力學(xué)學(xué)報(bào),2008,40(2)：255-260.

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

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

[11] Gerald C F,Wheatley PO.應(yīng)用數(shù)值分析[M].北京：機(jī)械工業(yè)出版社,2006.

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

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