榮俊杰,尤軍峰,文立華,校金友

(1.西北工業(yè)大學航天學院, 710072, 西安; 2.中國航天科技集團公司四院四十一所固體火箭發(fā)動機燃燒、熱結(jié)構(gòu)與內(nèi)流場國防科技重點實驗室, 710025, 西安)

彈性動力學瞬態(tài)問題的頻域邊界元快速解法

榮俊杰1,尤軍峰2,文立華1,校金友1

(1.西北工業(yè)大學航天學院, 710072, 西安; 2.中國航天科技集團公司四院四十一所固體火箭發(fā)動機燃燒、熱結(jié)構(gòu)與內(nèi)流場國防科技重點實驗室, 710025, 西安)

為解決常規(guī)基于離散傅里葉變換的頻域邊界元法難以解決無阻尼和低阻尼系統(tǒng)瞬態(tài)分析的問題,將指數(shù)窗口法與頻域邊界元法相結(jié)合,并采用預(yù)校正快速傅里葉變換(pFFT)方法加速邊界元求解。為進一步提高分析效率,針對頻域邊界元法所形成的系列線性方程組,提出了一種最小二乘外推法以獲得較高精度的迭代初值,可使初始解殘差小于10-2,從而顯著減少了迭代次數(shù);將新型的子空間回收算法用于頻域系列線性方程組的求解,加快了方程組的迭代收斂速度。算例表明,所提出的方法可顯著減少頻域邊界元法的迭代次數(shù),從而提高了計算效率,并有效降低了迭代解法的內(nèi)存消耗。

彈性動力學;邊界元法;迭代解法;子空間回收

彈性動力學的波動方程是研究固體中機械波傳播的理論基礎(chǔ),它的數(shù)值求解是力學、材料學、巖土工程學、地震學等眾多學科和工程領(lǐng)域的關(guān)鍵問題。邊界元方法是一種重要的工程數(shù)值方法,它在彈性動力學波動方程的數(shù)值求解中也得到了廣泛的應(yīng)用[1]。

彈性動力學問題的邊界元分析有3種基本途徑[2]:一是時域邊界元法,它采用時域基本解,同時在時間域和空間域進行離散求解;二是邊界-區(qū)域積分方程法,它采用差分法離散時間導(dǎo)數(shù)項,再利用彈性靜力學基本解獲得原問題對應(yīng)的邊界-區(qū)域積分方程;三是頻域邊界元法,它通過拉普拉斯或傅里葉變換將時域問題轉(zhuǎn)化為頻域中的橢圓邊值問題,采用邊界元法來求解,再通過逆變換獲得時域解。前2種方法存在時域積分的穩(wěn)定性問題,目前仍是國際研究的熱點。最近,時域卷積求積法被用于時域邊界元分析中,可以改善時間積分穩(wěn)定性問題,但其求解效率仍需進一步研究改進[3]。邊界-區(qū)域積分方程法需要計算體積分,可以采用對偶互易法、徑向積分法等,但計算量遠比邊界積分的大。

頻域邊界元法是一種較為簡便的動力學邊界元分析方法,它將時域問題轉(zhuǎn)化為若干個相互獨立的頻域問題來求解,編程實現(xiàn)簡便,適合于并行計算,并可以同時滿足動力學系統(tǒng)頻域和時域分析的需要,因此受到研究者的廣泛關(guān)注[4]。但是,當頻域邊界元法用于無阻尼或低阻尼系統(tǒng)的分析時,時域響應(yīng)的誤差很大甚至得不到正確的結(jié)果,這是由離散傅里葉變換的周期性引起的。文獻[5]采用指數(shù)窗口法(EWM)有效地解決了上述問題,并采用預(yù)校正快速傅里葉變換(pre-corrected FFT,簡稱pFFT)方法對頻域邊界元法進行了加速,大大提高了彈性動力學邊界元分析的效率。

本文旨在進一步提高頻域邊界元法求解彈性動力學波動問題的計算效率。在頻域邊界元法中,對于每個采樣頻率ωk求解一個線性方程組A(ωk)·x(ωk)=b(ωk),目前這些方程組是采用迭代方法獨立求解的。由于頻域邊界元法形成的系數(shù)矩陣A(ω)和右端向量b(ω)都是頻率ω的光滑函數(shù),所以解向量x(ω)可以用關(guān)于ω的有理函數(shù)較好地逼近。據(jù)此,在解得xk-1,…,xk-M,xk=x(ωk)之后,可以利用這些結(jié)果通過插值獲得對xk的迭代初值猜測,以減少求解迭代次數(shù),節(jié)省求解時間。GCRO-DR[6]是最近提出的一種求解系列線性方程組的算法,它利用子空間回收的思想重用上個線性方程組的求解信息,因此特別適用于求解系數(shù)矩陣光滑變化的序列方程組。本文在文獻[5]的基礎(chǔ)上,提出了解向量x(ω)的有理插值函數(shù)外插方法,并將GCRO-DR算法用于頻域邊界元線性方程組的求解,算例證明本文方法可以顯著減少迭代次數(shù)。

1 改進的動力學頻域邊界元法

設(shè)彈性區(qū)域為Ω,其邊界為Γ。彈性動力學問題的頻域邊界積分方程為

(1)

式中:ω為角頻率;cij為自由項;uj和σj分別為邊界位移和力;Uij和Tij分別為頻域彈性動力學基本解;左邊的積分為Cauchy主值積分。設(shè)畸變波和膨脹波速度分別為cp和cs,波數(shù)為kp和ks。Uij和Tij的表達式以及與積分方程(1)相對應(yīng)的邊界條件可參見文獻[2]。

將邊界Γ劃分為Ne個三角形單元,采用常元離散邊界積分方程(1)并考慮邊界條件,可得線性方程組

A(ω)x(ω)=b(ω)

(2)

式中:A(ω)為N×N的系數(shù)矩陣,N=3Ne;b(ω)為已知的N維列向量;x(ω)為包含未知位移和邊界力的待求向量。

為解決上述問題,文獻[7]將數(shù)字信號處理領(lǐng)域中的指數(shù)窗口法(EWM)引入到結(jié)構(gòu)動力學分析中。文獻[5]將該方法與頻域邊界元法相結(jié)合,有效地解決了常規(guī)頻域邊界元法存在的上述問題。

(3)

(4)

圖1 指數(shù)窗口法示意圖

設(shè)頻率和時間分辨率分別為Δω和Δt,則有

(5)

改進的頻域邊界元法的基本求解過程如下。

(6)

式中:ωk=kΔω-iη;k=0,1,…,Nω-1。

(2)分別對前(Nω/2+1)個頻率ωk(k=0,1,…,Nω/2)求解邊界積分方程(1),獲得頻域響應(yīng)xew(ωk),并由頻域響應(yīng)的共軛對稱性獲得后(Nω/2-1)個頻率上的響應(yīng)

式中:xew(k)=xew(kΔω-iη);x*表示x的復(fù)共軛。

(7)

(8)

2 線性方程組的加速求解

頻域邊界元法的主要計算量是對采樣頻率ωk=kΔω-iη求解系列方程組(2),或?qū)懗?/p>

Akx=bk,k=0,1,…,Nω/2

(9)

式中:Ak=A(ωk);bk=b(ωk)。本文采用pFFT快速邊界元法進行頻域求解。下面將首先闡述快速邊界元法與方程組迭代解法相結(jié)合的思想;然后,為進一步降低求解系列方程組(9)的計算量,提出一種通過最小二乘外推得到迭代初值的方法;最后,采用一種新型系列方程組解法來加速迭代收斂速度。

2.1 迭代解法及pFFT加速方法

當方程組(9)的階數(shù)N較大時,通常采用迭代解法求解(詳見2.3節(jié))。迭代解法的核心是計算矩陣-向量積Ax。直接計算Ax的計算量與N2成正比,即計算量為O(N2),不適合大規(guī)模計算。近20多年來,人們采用快速邊界元法進行Ax的加速計算,取得了巨大成功?？焖龠吔缭ò凑諉卧g的距離,將系數(shù)矩陣A分解為近場和遠場部分

A=Anear+Afar

(10)

2.2 外推法產(chǎn)生初始解

由式(2)可知,解向量x是頻率ω的函數(shù)。可進一步假設(shè)x是ω的光滑函數(shù)。一般地,x關(guān)于ω并不是處處光滑的,但多數(shù)采樣點附近的x仍是光滑的,因此該假設(shè)對多數(shù)采樣頻率ωk仍是成立的,這將在第3節(jié)通過數(shù)值算例予以證明。

由于x是ω的光滑函數(shù),當ω變化不大(即Δω較小)時,臨近的若干解xk=x(ωk)幾乎是線性相關(guān)的。于是,第k個方程組的解xk可以由前面K個解xk-i(i=1,2,…,K)的線性組合來獲得,即

(11)

式中:yi為待定系數(shù),可以通過將式(11)代入方程組Akxk=bk,由最小二乘法求得。令X為以前面K個解向量為列的矩陣

并令y=(y1,y2,…,yK)T,則式(11)可寫成

xk=Xy

(12)

分別用xk-i(i=1,2,…,K)和矩陣Ak相乘,可得矩陣S=AkX。于是,y是最小二乘問題

的解。采用QR分解方法,令QR=S,則y=R-1Qbk。于是

xk=XR-1Qbk

(13)

實際上,由于Δω不滿足足夠小的條件,而且x不一定是光滑函數(shù),按式(13)得到的解xk可能存在較大誤差,但卻可以作為一個較好的迭代初值,即令

(14)

上述產(chǎn)生迭代初值的方法與所采用的迭代解法無關(guān)。該方法在電磁場分析領(lǐng)域已有應(yīng)用,但在邊界元動力學分析中還未見采用。已有算例證明,用該方法獲得的初始解殘差在10-2以下,可以顯著減少總迭代次數(shù)。

2.3 系列線性方程組的高效解法

由以上基于系列方程組(9)的解向量之間的相關(guān)性,通過外推可以獲得較好的初始解,從而有效減少迭代次數(shù)。實際上,方程組(9)的系數(shù)矩陣Ak(k=0,1,…,Nω/2)之間也存在某種相關(guān)性(或相似性)。本文采用的GCRO-DR算法[6]是一種基于子空間回收思想的系列線性方程組解法,它利用系列方程組中相鄰系數(shù)矩陣Ak之間的相似性,每求解一個方程組,就從當前搜索空間中選擇一個近似不變子空間作為“回收”子空間,將其繼承到下一個方程組的解空間中?；厥兆涌臻g中包含了與收斂相關(guān)的信息,可以加快下一個方程組的迭代收斂速度。

設(shè)m為解的搜索空間最大維數(shù)。解第k個方程時,GCRO-DR首先在回收空間range(Us)(s

在頻域邊界元法中,隨著頻率的增大,解方程組(9)所需的迭代次數(shù)也會增加。在GCRO-DR算法中,從第(i-1)個方程繼承的回收子空間Us在用于求解第i個方程時,要計算AiUs,導(dǎo)致額外的s次矩陣向量乘法運算,限制了回收子空間的維數(shù)s。尤其當頻率較小、總迭代次數(shù)不大時,AiUs的額外計算量可能抵消子空間回收的優(yōu)勢。如果選擇較小的s,則在高頻時不能回收足夠的信息,會影響子空間回收的效率。為此,本文提出根據(jù)總迭代次數(shù)來動態(tài)調(diào)整s的大小,即選取s為當前方程組總迭代次數(shù)的β倍。大量數(shù)值算例表明,β的取值大致在0.15～0.25之間時,子空間回收的效率較高。s過大會使內(nèi)存和每步迭代的計算量增大,因此應(yīng)限定其上限sm。

3 算 例

采用文獻[5]中的算例1和算例2來驗證本文方法的計算精度和效率。分別采用3種線性方程組迭代解法:①全局GMRES算法[8](使用外推法求初值);②GCRO-DR算法(初值為零向量);③本文算法(采用GCRO-DR算法,用外推法求初值)。

將GCRO-DR算法和本文算法進行對比,可體現(xiàn)外推法的效率。全局GMRES算法是求解單個線性方程組時收斂速度最快的算法,將它與GCRO-DR算法進行比較,可進一步驗證子空間回收的優(yōu)勢。所有計算中均使用對角塊矩陣對系數(shù)矩陣進行預(yù)處理。

算例1階躍載荷下的棱柱

如圖2所示,棱柱長度L=3 m,左端固定,用階躍載荷p(t)加載,p0=-1 MN/m2。令泊松比為0,彈性模量E=211 GPa,密度ρ=7.85Mg/m3。計算中,將柱體表面劃分為3 200個單元。采樣周期T=0.0155s,采樣點數(shù)Nω=120,κ=3.0。GCRO-DR算法的參數(shù)m=70,sm=20,外推法的參數(shù)K=4,收斂誤差為1×10-7。

圖2 階躍載荷下的棱柱棒

首先考察本文算法的精度。在設(shè)定泊松比為0的條件下,此問題具有解析解[5]。圖3給出了固定端應(yīng)力隨時間的變化關(guān)系。與解析解對比,本文方法得到的數(shù)值解與之基本一致。

t:實驗時間; c:波速; L:棱柱長度

為驗證本文提出的參數(shù)k選取方法的合理性,分別取s=5,25和β=0.2,使用本文算法進行計算,得到迭代次數(shù)隨頻率的變化曲線,如圖4所示。當s取5時,在高頻范圍內(nèi)迭代次數(shù)較多,而當s取20時,在低頻范圍內(nèi)迭代次數(shù)較多。只有當β=0.2時,在全部頻域范圍內(nèi)迭代次數(shù)都基本處于最低水平。

圖4 本文算法的迭代次數(shù)隨頻率的變化曲線

為驗證本文算法的效率,分別用上述3種算法對算例進行計算,表1給出了3種算法的總迭代次數(shù)。與全局GMRES算法相比,本文算法的迭代次數(shù)有所減少,但減少量有限,這是由于本文算法所使用的GCRO-DR算法在本質(zhì)上是對重啟型GMRES的改進。隨著迭代次數(shù)的增加,全局GMRES算法的內(nèi)存消耗及每步迭代的計算量都會增大。表1列出的內(nèi)存使用情況表明,本文算法的內(nèi)存消耗小于全局GMRES算法的消耗。與GCRO-DR算法相比,本文算法由于使用了本文提出的外推法計算迭代初值,明顯減少了迭代次數(shù)。

表1 3種算法計算算例1的總迭代次數(shù)和 內(nèi)存占用量

算法求解所有方程的總迭代次數(shù)占用內(nèi)存/MB全局GMRES算法551524 8GCRO?DR算法648717 8本文算法480218 6

算例2沖擊載荷作用下的帶孔平板

圖5 沖擊載荷下的帶孔平板

通過一個更實際的問題來驗證本文算法的精度和效率。圖5給出了一鋁制平板結(jié)構(gòu)的尺寸、載荷分布以及約束情況。鋁板的材料參數(shù)分別為:E=69 GPa,ρ=2.7 Mg/m3,ν=0.3。假定沖擊載荷均勻施加在平板的頂部,平板的底部固定,如圖5所示。試樣表面被劃分為14 072個三角單元。頻域邊界元法的參數(shù)為:T=0.003,Nω=128,κ=2.0。分別使用上述3種算法進行計算。GCRO-DR算法的參數(shù)選擇為m=40,sm=20,β=0.25,外推法參數(shù)K=4,收斂誤差為1×10-4。對于圖5中S點的應(yīng)變,本文算法的計算結(jié)果如圖6所示,可見除去一些不確定因素外,本文算法所得的結(jié)果與試驗結(jié)果相符。表2給出了3種算法的總迭代次數(shù)和內(nèi)存占用量,從中可以看出,本文算法與其他算法相比迭代次數(shù)顯著減少,與全局GMRES算法相比,總迭代次數(shù)減少了57%,占用內(nèi)存減少了56%。

圖6 S點應(yīng)變隨時間變化圖

表2 3種算法計算算例2的總迭代次數(shù)和 內(nèi)存占用量

算法求解所有方程的總迭代次數(shù)占用內(nèi)存/MB全局GMRES算法8122132 4GCRO?DR算法659657 4本文算法350861 3

4 結(jié) 論

本文將指數(shù)窗口法與頻域邊界元法相結(jié)合,有效地解決了常規(guī)頻域邊界元法存在的難以求解無阻尼、低阻尼系統(tǒng)以及求解效率低的問題。頻域邊界元分析采用pFFT方法進行加速,實現(xiàn)簡便,加速效率高。為了進一步提高頻域邊界元法的計算效率,針對頻域邊界元法所形成的系列線性方程組,本文提出了一種最小二乘外推法,用于獲得較高精度的迭代初值,同時將文獻[6]提出的子空間回收算法用于頻域系列方程組的快速求解,并給出了其控制參數(shù)的取值方法。數(shù)值算例表明,本文算法可以顯著降低頻域邊界元分析的迭代次數(shù),并有效降低迭代解法的內(nèi)存占用量。

[1] 姚振漢, 王海濤.邊界元法 [M].北京: 高等教育出版社, 2010: 162-189.

[2] 稽醒, 臧躍龍, 程玉民.邊界元法進展及通用程序 [M].上海: 同濟大學出版社, 1997: 51-101.

[3] BANJAI L, SAUTER S.Rapid solution of the wave equation in unbounded domains [J].SIAM Journal on Numerical Analysis, 2008, 47(1): 227-249.

[4] CHAILLAT S, BONNET M, SEMBLAT J F.A multi-level fast multipole BEM for 3-D elastodynamics in the frequency domain [J].Computer Methods in Applied Mechanics and Engineering, 2008, 197(49): 4233-4249.

[5] XIAO Jinyou, YE Wenjing, CAI Yaxiong, et al.Precorrected FFT accelerated BEM for large-scale transient elastodynamic analysis using frequency-domain approach [J].International Journal for Numerical Methods in Engineering, 2012, 90(1): 116-134.

[6] PARKS M L, STURLER E, MACKEY G, et al.Recycling Krylov subspaces for sequences of linear systems [J].SIAM Journal on Scientific Computing, 2006, 28(5): 1651-1674.

[7] KAUSEL E, ROESSET J M.Frequency domain analysis of undamped systems [J].Journal of Engineering Mechanics, 1992, 118(4): 721-734.

[8] SAAD Y, MARTIN H S.GMRES: a generalized minimal residual algorithm for solving nonsymmetric linear systems [J].SIAM Journal on Scientific Computing, 1986, 7(3): 856-869.

[本刊相關(guān)文獻鏈接]

李應(yīng)剛,陳天寧,王小鵬,等.外部動態(tài)激勵作用下齒輪系統(tǒng)非線性動力學特性.2014,48(1):101-105.[doi:10.7652/xjtuxb201401017]

周生喜,曹軍義,Alper ERTURK,等.壓電磁耦合振動能量俘獲系統(tǒng)的非線性模型研究.2014,48(1):106-111.[doi:10.7652/xjtuxb201401018]

張靜,郭宏偉,劉榮強,等.空間含鉸可展桁架結(jié)構(gòu)的非線性動力學建模與分析.2013,47(11):113-119.[doi:10.7652/xjtuxb201311020]

劉俊俏,苗福生,李星.二維各向異性功能梯度材料熱傳導(dǎo)的邊界元分析.2013,47(5):77-81.[doi:10.7652/xjtuxb2013 05014]

王燾,校金友,曹衍闖,等.采用Galerkin離散方法的T-小波邊界元法.2010,44(12):99-104.[doi:10.7652/xjtuxb2010 12019]

(編輯 葛趙青)

AFastFrequencyDomainBoundaryElementMethodforTransientElastodynamicAnalysis

RONG Junjie1,YOU Junfeng2,WEN Lihua1,XIAO Jinyou1

(1.College of Astronautics, Northwestern Polytechnical University, Xi’an 710072, China; 2.National Key Laboratory of Combustion, Flow and Thermo-Structure, The 41st Institute of the Forth Academy of CASC, Xi’an 710025, China)

The conventional frequency-domain boundary element method (BEM) based on discrete Fourier transform is difficult to complete transient analysis of low damped and undamped systems.To solve this problem, the exponential window method is combined with the frequency-domain BEM, and the precorrected fast Fourier transform (pFFT) is used to accelerate the frequency-domain BEM analysis.For solving the sequence of linear systems corresponding to the sampling frequencies with higher efficiency, an extrapolation method is proposed to obtain good initial guesses.In addition, a newly developed Krylov subspace recycling method is employed to solve this sequence of linear systems.Numerical examples demonstrate a dramatic reduction of iterations, which leads to higher efficiency and smaller memory consuming.

elastodynamics; boundary element method; iterative method; Krylov subspace recycling

10.7652/xjtuxb201403021

2013-05-07。

榮俊杰(1990—),男,碩士生;校金友(通信作者),男,副教授。

國家自然科學基金資助項目(11102154, 11074201);教育部高等學校博士學科點專項科研基金資助項目(2010610212009, 2011610211006)。

時間: 2013-12-25

O353.2

:A

:0253-987X(2014)03-0115-06

網(wǎng)絡(luò)出版地址: http:∥www.cnki.net/kcms/detail/61.1069.T.20131225.1701.004.html