王博,王燾,尤軍峰,校金友

(1.西北工業(yè)大學(xué)航天學(xué)院,710072,西安;2.中國航天科技集團公司四院四十一所,710025,西安)

復(fù)雜薄壁結(jié)構(gòu)的聲振耦合分析屬于多物理場耦合的科學(xué)計算問題,也是諸多工業(yè)部門對其產(chǎn)品進行力學(xué)環(huán)境預(yù)示的核心問題,具有非常重要的科學(xué)意義和極其廣泛的工程應(yīng)用價值。尤其是在航天領(lǐng)域,航天器多為薄壁結(jié)構(gòu),剛度小,而在發(fā)射過程中產(chǎn)生的噴氣噪聲和氣動噪聲強度大、頻帶寬,通?？蛇_10～10 000 Hz,再考慮到火箭發(fā)動機燃燒不均勻產(chǎn)生的振動環(huán)境,聲振耦合現(xiàn)象尤為明顯[1]。

簡單結(jié)構(gòu)的聲振耦合分析可以通過解析方法得到[2],但對于復(fù)雜結(jié)構(gòu)則必須借助數(shù)值手段[3-4]。聲振耦合問題的數(shù)值解法主要包括有限元法(FEM)、邊界元-有限元(BEM-FEM)耦合方法和統(tǒng)計能量法(SEA)。其中,BEM-FEM耦合方法能夠充分發(fā)揮BEM和FEM的各自優(yōu)勢,即對聲場和彈性結(jié)構(gòu)分別采用BEM和FEM進行建模,再根據(jù)結(jié)構(gòu)與聲場界面上的邊界條件,將聲學(xué)BEM方程和結(jié)構(gòu)FEM方程聯(lián)立求解,得到聲場和結(jié)構(gòu)的響應(yīng),是對聲振耦合問題進行仿真計算的理想工具。

構(gòu)造基于BEM-FEM的聲振耦合方程需要完成聲場網(wǎng)格和結(jié)構(gòu)網(wǎng)格之間的數(shù)據(jù)交換。匹配網(wǎng)格間的數(shù)據(jù)交換可以直接通過界面上相互作用的強形式來準(zhǔn)確獲得[5]。由于聲場和應(yīng)力場的場函數(shù)梯度不同,采用不同的單元尺寸和形函數(shù)來離散結(jié)構(gòu)和聲場通常能獲得更高的分析效率,同時也能夠兼顧建模的便利性。對于聲振耦合問題,數(shù)據(jù)交換的準(zhǔn)確程度直接影響到系統(tǒng)方程的求解精度和收斂性[6]。目前,非匹配網(wǎng)格間的數(shù)據(jù)交換算法主要分為插值類[7]、徑向基函數(shù)類[8]和加權(quán)殘差類[9]等,其中加權(quán)殘差類方法使非匹配耦合界面上的邊界條件以弱形式成立,考慮了虛功等效原理,能夠確保傳遞過程的能量守恒[7],當(dāng)使用高階單元時,比其他2種算法具有更高的精度[10]。

雖然BEM-FEM耦合方法是研究聲振耦合問題的一種有效手段,但傳統(tǒng)邊界元法形成的系數(shù)矩陣是滿陣,其計算復(fù)雜度為O(N2)(N為系統(tǒng)的自由度數(shù)目),難以直接用于分析大規(guī)模問題。近30年來發(fā)展起來的快速算法,如快速多極方法、H矩陣法、自適應(yīng)交叉近似方法和小波邊界元法等[11],利用系數(shù)矩陣可分塊低秩分解的特性,將計算量降低到了O(NlogN),推動了BEM在實際工程領(lǐng)域的應(yīng)用。

BEM快速算法的引入為通過BEM-FEM耦合方法分析大規(guī)模聲振耦合問題提供了有力支撐,然而目前相關(guān)方面的研究大多集中在低頻段[4]。這是因為,對于高頻波動問題,采用上述這些快速算法得到的系數(shù)矩陣子矩陣將不再能夠進行低秩分解[12],因此這些方法只能用于加速計算低頻問題。2006年提出的寬頻快速多級算法[12]采用將快速多級算法中變換矩陣對角化的方法,將高頻波動問題的計算量降低到了O(NlogN),但是這種方法的實現(xiàn)過程十分復(fù)雜。近年來基于高頻問題定向低秩特性提出的快速定向算法[13],同樣成功地將計算量降低到了O(NlogN)。相比于寬頻快速多級算法,該算法不需要對積分核函數(shù)進行復(fù)雜的解析展開,構(gòu)造低秩分解的方法靈活,實現(xiàn)過程簡單,而且與核函數(shù)的具體表達式無關(guān),因而得到了廣泛關(guān)注和迅速發(fā)展。但是,目前該方法僅用于處理單一的聲學(xué)問題,在聲振耦合問題分析中的應(yīng)用還未見報道。

本文首先采用Burton-Miller方程描述聲場,并通過二階Nystr?m方法對聲場進行離散;接著,考慮更具一般意義的耦合界面上聲場網(wǎng)格和結(jié)構(gòu)網(wǎng)格不匹配的情形,同時結(jié)合本文所選聲場和結(jié)構(gòu)網(wǎng)格的特點,基于加權(quán)余量法建立耦合界面上聲場和結(jié)構(gòu)網(wǎng)格之間的數(shù)據(jù)交換矩陣;然后,采用雙向耦合方法將聲場BEM方程與結(jié)構(gòu)FEM方程聯(lián)立,建立聲振耦合系統(tǒng)方程,并給出其求解格式;最后,通過適用于寬頻問題的快速定向算法對耦合方程中的常規(guī)BEM矩陣進行加速,實現(xiàn)了對大規(guī)模、寬頻聲振耦合問題的快速分析,并通過三維數(shù)值算例驗證了本文方法的正確性以及相較現(xiàn)有同類方法的優(yōu)越性。

1 BEM-FEM耦合方法

為了獲得更高的分析精度,采用雙向耦合方法對聲振耦合系統(tǒng)進行建模,即將耦合界面上聲場的聲壓視為附加節(jié)點載荷作用于結(jié)構(gòu),同時將耦合界面上的結(jié)構(gòu)振動響應(yīng)視為附加聲激勵作用于聲場,最終獲得雙向耦合方程。為簡便起見,本節(jié)僅建立聲場網(wǎng)格和結(jié)構(gòu)網(wǎng)格完全匹配情形下的聲振耦合方程,在第2節(jié)再對非匹配情形作進一步討論。

1.1 建立BEM-FEM耦合方程

對結(jié)構(gòu)采用FEM建模,考慮了耦合界面上聲壓作用的頻域結(jié)構(gòu)動力學(xué)方程如下

(K-ω2M+iωC)u=fs+Ta2sp

(1)

式中:K、M和C分別為剛度矩陣、質(zhì)量矩陣和阻尼矩陣;u為結(jié)構(gòu)FEM網(wǎng)格上節(jié)點的位移向量;fs為節(jié)點的載荷向量;p為聲場BEM網(wǎng)格高斯點上的聲壓向量;Ta2s為耦合矩陣,作用是將p轉(zhuǎn)化成結(jié)構(gòu)網(wǎng)格節(jié)點上的等效節(jié)點力,由于聲場網(wǎng)格和結(jié)構(gòu)網(wǎng)格完全匹配,故Ta2s僅相當(dāng)于一個對面內(nèi)分布載荷p的等效節(jié)點力的轉(zhuǎn)化矩陣,通過單元積分即可獲得。

采用二階Nystr?m方法離散聲學(xué)Burton-Miller方程,得到聲學(xué)BEM方程如下

Hp=Gq

(2)

式中:H、G為BEM系數(shù)矩陣;q為聲壓法向?qū)?shù)向量。

根據(jù)線性聲學(xué)理論,聲壓p與聲場速度勢Φ存在如下關(guān)系

p=iρωΦ

(3)

式中ρ為聲介質(zhì)密度。對式(3)沿BEM網(wǎng)格外法向na求導(dǎo),可得

q=iρωva·na=ρω2ua·na

(4)

式中:va為聲場網(wǎng)格高斯積分點的法向振速向量;ua為相應(yīng)的位移向量。由式(4)可知,結(jié)構(gòu)網(wǎng)格表面的位移u經(jīng)變換可以間接表示式(2)中的q,即

q=Ts2au

(5)

式中Ts2a為耦合矩陣。

將式(5)代入式(2),經(jīng)整理可得

-GTs2au+Hp=0

(6)

定義結(jié)構(gòu)動剛度矩陣Ks為

Ks=K-ω2M+iωC

(7)

將式(7)代入式(1),整理后可得

Ksu-Ta2sp=fs

(8)

聯(lián)立式(6)和式(8)并寫成矩陣形式,可得

(9)

式(9)即為基于BEM-FEM的聲振耦合系統(tǒng)方程。

1.2 耦合線性方程組求解

對于式(9)描述的BEM-FEM耦合方程,其系數(shù)矩陣同時包含了FEM矩陣、BEM矩陣和數(shù)據(jù)交換矩陣,各子塊矩陣的性質(zhì)也完全不同。一般而言,FEM動剛度矩陣Ks和耦合矩陣Ta2s為稀疏矩陣,而BEM矩陣G和H為稠密矩陣,見圖1。

圖1 基于BEM-FEM的聲振耦合系統(tǒng)矩陣非零元素分布

針對此類線性方程組,采用單一的求解方法通常難以獲得最優(yōu)的求解效率,有時甚至難以解出可靠結(jié)果,因此必須針對該系數(shù)矩陣的特點,設(shè)計專門的求解格式,并為系數(shù)矩陣的不同子塊選擇最為合適的求解器,以便快速獲得高精度的計算結(jié)果。

采用耦合的Burton-Miller方法求解式(9),將結(jié)構(gòu)系統(tǒng)FEM動剛度矩陣聚縮進聲學(xué)BEM系統(tǒng)矩陣中,從而得到耦合的Burton-Miller方程,再迭代求解該方程得到聲壓,最后將得到的聲壓代入結(jié)構(gòu)FEM方程得到位移解。這種方法首先通過消去變量,減少了方程組維數(shù),具有較高的計算效率[4]。

對式(8)作如下變換

(10)

將式(10)代入式(6),經(jīng)整理即可得到如下耦合的Burton-Miller方程

(11)

求解式(11)得到聲壓p,再將其代入式(10)計算得到位移u,最終獲得整個聲振耦合系統(tǒng)的響應(yīng)。

分析上述過程可以發(fā)現(xiàn),式(11)的求解效率和精度決定了整個聲振耦合系統(tǒng)的分析效率和精度。

在式(11)中,FEM方程的求解屬于稀疏線性方程組求解問題,應(yīng)該采用直接解法。本文采用多波前LU分解算法求解FEM方程,通過引入一些典型的修正策略,如特殊的選主元方法、對系數(shù)矩陣元素的重排序以及巧妙地利用數(shù)據(jù)存儲結(jié)構(gòu)等,可使該方法實現(xiàn)近似線性的計算復(fù)雜度。

為了加速GMRES算法的迭代收斂過程,對系數(shù)矩陣進行合適的預(yù)處理十分必要。文獻[4]已經(jīng)證明,當(dāng)采用直接解法求解式(9)中的FEM方程時,對整個耦合矩陣進行預(yù)處理相當(dāng)于僅對邊界元矩陣H進行預(yù)處理。常見的BEM矩陣預(yù)處理方法主要有塊對角方法和ILU方法,前者實現(xiàn)簡單,內(nèi)存占用少,而后者通常可以獲得更高的整體效率。為了盡可能地快速求解式(9),本文選擇ILU方法。采用預(yù)處理技術(shù)處理后,最終獲得的耦合的Burton-Miller方程如下

(12)

2 非匹配網(wǎng)格數(shù)據(jù)交換

前面介紹了聲場網(wǎng)格與結(jié)構(gòu)網(wǎng)格匹配情形下的快速BEM-FEM聲振耦合系統(tǒng)方程構(gòu)造方法,但在實際應(yīng)用中,這種建模方法存在明顯的局限性:首先,不同物理場的場函數(shù)通常對單元尺寸有不同的要求,采用匹配網(wǎng)格離散必須考慮對單元尺寸要求最苛刻的情況,從而增加不必要的計算負擔(dān);其次,聲場和結(jié)構(gòu)場通常在不同的軟件環(huán)境中離散,確保網(wǎng)格的完全匹配比較困難,不利于建模。

基于此,本節(jié)討論非匹配網(wǎng)格數(shù)據(jù)交換技術(shù)在快速BEM-FEM聲振耦合方程中的應(yīng)用,并針對所采用單元的特點,利用加權(quán)余量法發(fā)展了一種非匹配高階協(xié)調(diào)元和高階非協(xié)調(diào)元之間的數(shù)據(jù)交換算法,該算法能確保傳遞過程的能量守恒,因此,相較插值類方法和徑向基函數(shù)類方法可以獲得更高的數(shù)據(jù)交換精度。

將聲場邊界離散為二次曲面6節(jié)點三角形單元,其單元積分采用6點高斯積分公式;將結(jié)構(gòu)離散為二次曲面8節(jié)點四邊形單元,其單元積分采用9點高斯積分公式。單元示意圖見圖2。

(a)結(jié)構(gòu)FEM單元(b)聲場BEM單元●:節(jié)點; +:高斯積分點圖2 單元示意圖

首先,建立能量守恒方程,此時需要在聲場Ψa和結(jié)構(gòu)場Ψs之間引入一個虛擬物理場Ψv,并分別在Ψv和Ψs之間、Ψv和Ψa之間定義如下形式的能量對Evs和Eva

(13)

(14)

式中:Nv、Ns和Na分別為虛擬物理場網(wǎng)格Ωv、結(jié)構(gòu)網(wǎng)格Ωs和聲場網(wǎng)格Ωa對應(yīng)的形函數(shù)矩陣;Vv、Vs和Va分別為定義在Ωv、Ωs和Ωa節(jié)點上的場向量。令Evs和Eva相等,即

(15)

定義提供數(shù)據(jù)的場為源場,接收數(shù)據(jù)的場為目標(biāo)場。當(dāng)結(jié)構(gòu)為源場、聲場為目標(biāo)場時,令Ωv與Ωa重合,從而得到Nv=Na,代入式(15)并經(jīng)整理可得

(16)

定義如下2個數(shù)據(jù)交換矩陣

(17)

將式(17)代入式(16)可得

(18)

式(18)即為結(jié)構(gòu)網(wǎng)格節(jié)點向量Vs向聲場網(wǎng)格節(jié)點向量Va轉(zhuǎn)化的公式,可以發(fā)現(xiàn)當(dāng)聲場網(wǎng)格和結(jié)構(gòu)網(wǎng)格完全匹配時,Va=Vs。同理,當(dāng)聲場為源場、結(jié)構(gòu)為目標(biāo)場時,聲場網(wǎng)格節(jié)點向量Va向結(jié)構(gòu)網(wǎng)格節(jié)點向量Vs轉(zhuǎn)化的公式為

(19)

式中,數(shù)據(jù)交換矩陣Dss和Dsa通過下式計算

(20)

接下來,分別推導(dǎo)耦合矩陣Ts2a和Ta2s的具體表達式。Ts2a的功能是將結(jié)構(gòu)網(wǎng)格節(jié)點上的位移轉(zhuǎn)化成聲場網(wǎng)格高斯積分點上的法向振速,Ta2s的功能是將聲場網(wǎng)格高斯積分點上的聲壓轉(zhuǎn)化成結(jié)構(gòu)網(wǎng)格節(jié)點上的等效節(jié)點力。

(21)

為了獲得聲場網(wǎng)格高斯積分點處的位移向量uaG,對uav進行單元插值

uaG=Na·uav

(22)

對式(22)中的uaG沿法向求導(dǎo),即可獲得聲場網(wǎng)格高斯積分點處的法向振速

qaG=Nna·uaG

(23)

式中,Nna為聲場網(wǎng)格高斯積分點處的法向求導(dǎo)矩陣。整理上述過程,可得

(24)

式(24)即為Ts2a的具體表達式。

由于本文采用非協(xié)調(diào)元離散聲場邊界,故p為聲場網(wǎng)格高斯積分點上的聲壓向量,進行數(shù)據(jù)交換前,首先需對p進行一次外插,以獲得聲場網(wǎng)格節(jié)點上的聲壓向量

(25)

(26)

為了將聲壓作為附加節(jié)點力載荷施加到結(jié)構(gòu)節(jié)點上,先對其進行單元插值,獲得結(jié)構(gòu)網(wǎng)格高斯積分點上的聲壓分布

psG=Ns·psv

(27)

然后采用高斯積分計算單元內(nèi)分布載荷的等效節(jié)點力轉(zhuǎn)換公式,得到結(jié)構(gòu)網(wǎng)格節(jié)點上的等效節(jié)點力

(28)

式中,Ws由高斯積分權(quán)重與相應(yīng)的方向?qū)?shù)的乘積組成。整理式(25)～(28),可得

(29)

式(29)即為Ta2s的具體表達式。

由上述推導(dǎo)過程可知,式(24)和(29)中對矩陣Daa和Dss進行求逆是這個方法最耗時的部分。不過,對式(17)和(20)進行分析可以發(fā)現(xiàn),這2個矩陣類似于FEM中的質(zhì)量矩陣,是對稱正定的,因此可以通過Cholesky分解方法實現(xiàn)快速分解,當(dāng)問題規(guī)模較大時,也可以借助一些快速Cholesky分解方法,如多波前Cholesky分解方法等。

式(17)和(20)定義了4個數(shù)據(jù)交換矩陣的計算格式,其中Daa和Dss的計算可以通過簡單的單元積分實現(xiàn)。但是,Das和Dsa的計算都涉及到對定義在不同網(wǎng)格上的形函數(shù)矩陣乘積進行積分的問題,需要對積分點進行投影,以及對單元重疊部分進行切割并進行重劃分等復(fù)雜操作,實現(xiàn)過程復(fù)雜,這也是非匹配網(wǎng)格數(shù)據(jù)交換中的重點和難點。

本文采用文獻[7]給出的基于單元的積分方案計算Das和Dsa。為了敘述簡便,后文稱被積單元為“主單元”,稱與主單元重合的單元為“從單元”。該方法只需確定出主單元積分點的投影點屬于哪個從單元,同時計算出投影點在從單元內(nèi)相應(yīng)的參數(shù)坐標(biāo)即可。遍歷全部主單元后,將前期確定的投影點形參坐標(biāo)代入式(17)和(20),對每一個主單元直接進行積分獲得其對應(yīng)的子矩陣,最后根據(jù)投影點所屬的從單元編號將子矩陣投放到Das和Dsa中。

采用基于單元的積分方案,主要工作量是確定主單元積分點的投影點所屬的從單元及其參數(shù)坐標(biāo)。投影點在從單元上的大體位置可以根據(jù)幾何模型的建模特點或快速BEM中的八叉樹狀結(jié)構(gòu)快速獲得,但具體形參坐標(biāo)的確定還存在一定困難,特別是針對高階單元,因此本文采用一種牛頓類迭代方法來確定投影點的形參坐標(biāo)。

3 快速定向邊界元法

通過常規(guī)BEM獲得的系數(shù)矩陣H和G均為滿陣,計算量和存儲量都很大,無法運用大規(guī)模問題的分析計算方法,因此本文采用快速定向壓縮邊界元法對其進行快速計算,可在寬頻范圍內(nèi)將計算復(fù)雜度降低到O(NlogN)。下面僅對這種方法作簡單介紹,詳細內(nèi)容可參考文獻[13]。

首先定義自適應(yīng)八叉樹結(jié)構(gòu):定義格子C0,使所有聲場邊界單元均包含在其中,作為八叉樹的根;然后對八叉樹中的葉格子進行細分,得到下一層的格子;繼續(xù)這樣細分,直到每一個葉格子都只包含不超過Np個聲場節(jié)點。對于每一層格子,將它們的尺寸記為w,將w<λ的層稱為低頻層,其他的層稱為高頻層。對于格子C,若它在低頻層上,則定義它的近場NC為同一層上所有與它相鄰的格子所在的區(qū)域;若它在高頻層上,則定義它的近場為所有格子{B|dist(B,C)≤O(kw2)}所在的區(qū)域。定義它的遠場FC為除去近場NC后的其他區(qū)域。將C的父格子記為P,定義它的交互場IC=NPNC。當(dāng)C在高頻層上時,將它的交互場劃分成若干個楔形區(qū),使得每個楔形區(qū)張開的角度為O(1/(kw)),則格子C和它的每個楔形區(qū)之間均滿足拋物線分離條件,如圖3所示。

圖3 拋物線分離條件

將圖3所示的楔形區(qū)稱為C的a方向楔形區(qū)。為便于敘述,將低頻層格子的遠場稱為它的(0,0,0)方向的楔形區(qū)?？梢宰C明:不論C在高頻層還是在低頻層,當(dāng)源點y和場點x分別在C中和它的楔形區(qū)中時,核函數(shù)G(x,y)均可按如下形式進行低階展開

(30)

因此,H、G中與C及其楔形區(qū)相對應(yīng)的子矩陣可以近似低秩分解,即與此子矩陣相對應(yīng)的乘法運算可以采用快速算法加速。

與傳統(tǒng)快速多極方法不同,本文采用的快速定向壓縮算法并不將核函數(shù)解析展開成式(30)的形式,即不對核函數(shù)進行多極展開或局部展開,而是借助于等效源和檢測勢的方法構(gòu)造快速算法?？焖俣ㄏ驂嚎s算法具有與傳統(tǒng)快速多極方法相同的結(jié)構(gòu),即都由S2M、M2M、M2L、L2L、L2T變換組成,下面對它們分別進行介紹。

S2M變換根據(jù)每個葉格子C中的源,構(gòu)造它的外出等效源。為此,首先定義外出等效點和外出檢測點,如圖4所示。S2M變換由如下2步運算構(gòu)成:

(1)根據(jù)C中節(jié)點上的源計算外出檢測點上的外出檢測勢,其運算矩陣記為Eup;

S2M運算過程如圖4所示。

若C為非葉格子,M2M運算根據(jù)它子格子上的外出等效源計算它本身的外出等效源。此運算步驟與S2M類似,只需要將第(1)步中的源換為子格子的等效源即可。需注意,當(dāng)C在高頻層上時,需要對它的每一個楔形區(qū)構(gòu)造一組等效源。

(a)低頻層上的S2M變換

(b)高頻層上的S2M變換圖4 S2M變換過程

(a)低頻層上的L2T變換

在快速定向壓縮算法的下行運算中,對于每一個格子需要定義另外2組點,稱為返回等效點和返回檢測點,它們可以定義為分別與外出檢測點和外出等效點重合,如圖5所示。類似地,快速算法中將計算返回檢測點上的勢,稱為返回檢測勢,并且需要在返回等效點上構(gòu)造等效源,稱為返回等效源。

(b)高頻層上的L2T變換圖5 L2T變換過程

快速定向壓縮算法的M2L運算,是根據(jù)格子C上的外出等效源計算IC中與C在同一層格子上的下行檢測勢。當(dāng)格子B∈IC且B與C在同一層時,若B在C的a方向楔形區(qū)中,則C必在B的-a方向楔形區(qū)中。因此,格子B上a方向的返回檢測勢需要用格子C上-a方向的外出等效源計算。

對于每一個非空、非葉格子C,L2L運算是根據(jù)它的返回檢測勢計算它子格子的返回檢測勢,而對于葉格子,L2T運算是根據(jù)它的返回檢測勢計算其內(nèi)部配置點上的勢。L2L和L2T可分別看作M2M和S2M的逆變換。L2T的計算過程如圖5所示。

4 算 例

全部算例均在一臺12核(Intel Xeon E5-2643 v4,3.40 GHz)128 GB內(nèi)存的工作站上進行,所有程序都采用單線程C++編寫。BEM快速算法的精度和GMRES算法的收斂殘差都給定為10-6。為了驗證本文所提之方法的正確性與性能,首先引入2個具有代表性的聲振耦合解析模型。

模型1聲振耦合內(nèi)問題。一單位立方體內(nèi)聲腔,聲介質(zhì)為水,其中一個表面簡支了一塊1 m×1 m的鋼板,厚度t=0.01 m,其余表面固支,在鋼板表面距離其左下角(0.2 m,0.3 m)的位置作用一單位簡諧集中力,見圖6。水的密度ρa=1 000 kg/m3,水中聲速va=1 481 m/s;鋼的彈性模量E=210 GPa,密度ρs=7 900 kg/m3,泊松比μ=0.3。將該模型的結(jié)構(gòu)部分離散為5 400個四邊形殼單元,聲場部分離散為2 580個三角形殼單元,BEM-FEM模型的總自由度為37 428。

文獻[2]給出了該模型激振點處位移和聲壓的解析解。

模型2聲振耦合外問題。一鋼質(zhì)球殼,半徑R=5 m,壁厚t=0.05 m,完全浸入水中,在其表面任一點處作用一單位簡諧集中力,見圖7。水和鋼的物理參數(shù)同模型1。將該模型的結(jié)構(gòu)部分離散為2 400個8節(jié)點四邊形殼單元,聲場部分離散為4 800個6節(jié)點三角形殼單元,BEM-FEM模型的總自由度為72 012。

該模型激振點處位移和聲壓的解析解可參見文獻[2]。

圖6 模型1的FEM模型圖7 模型2的 FEM模型

4.1 數(shù)據(jù)交換精度驗證

此算例旨在驗證第2節(jié)構(gòu)造的數(shù)據(jù)交換算法的精度。定義如下試函數(shù)

f(x,y,z)=x2+y2+z2+yz+zx+xy+1

(31)

式中,(x,y,z)為耦合界面上的網(wǎng)格節(jié)點坐標(biāo)。由于聲場和結(jié)構(gòu)均采用二階單元離散,故試函數(shù)取二次即可。

(32)

(a)ζs2a

(b)ζa2s圖8 模型1的數(shù)據(jù)交換精度驗證

(a)ζs2a

(b)ζa2s圖9 模型2的數(shù)據(jù)交換精度驗證

觀察圖8和圖9可知:當(dāng)模型外表面的幾何形狀為平面時,數(shù)據(jù)交換過程不損失精度;當(dāng)外表面為曲面時,部分節(jié)點的最大誤差為10-3量級,但整體誤差都在10-4以下。分析圖9中誤差產(chǎn)生的原因,認為主要有2個:①對曲面進行積分時受本文所采用積分方案的精度限制;②2套網(wǎng)格之間存在間隙,模型有離散誤差。

4.2 聲振耦合精度驗證

此算例旨在驗證本文構(gòu)造的快速BEM-FEM聲振耦合方法的正確性及性能。分別對模型1和模型2進行掃頻分析,模型1的掃頻區(qū)間為1～1 000 Hz,步長為1 Hz;模型2的掃頻區(qū)間為1～300 Hz,步長取為1 Hz。分別提取這2個模型每個采樣頻率下激振點的位移,繪制頻響曲線,并和解析解進行對比,見圖10和圖11。

圖10 模型1激振點位移頻響曲線的數(shù)值解與解析解對比

圖11 模型2激振點位移頻響曲線的數(shù)值解與解析解對比

觀察圖10和圖11可以發(fā)現(xiàn),本文方法計算出的這2個模型的激振點位移均與解析解吻合良好,證明無論是對于聲振耦合內(nèi)問題還是外問題,本文方法均能保證正確性。在公開文獻中,模型1的掃頻范圍僅為1～250 Hz[14],而模型2的掃頻范圍僅為1～100 Hz[4,6],遠低于本文方法所適用的頻段。此外,將上述本文計算結(jié)果與相關(guān)文獻中同類方法的計算結(jié)果[4,14]進行對比,發(fā)現(xiàn)本文方法的精度也具有一定優(yōu)勢。

在計算效率方面,以模型2浸水球殼算例進行對比,計算頻率為65 Hz。文獻[4]采用同類方法,在Intel Xeon 5160 3.0 GHz GPU、單線程C++代碼的計算環(huán)境下計算一個總自由度為27 958(其中FEM自由度為23 964,BEM自由度為3 994)的模型,最少耗時3.46 min,而本文方法在類似計算環(huán)境下計算一個總自由度為72 012(其中FEM自由度為43 212,BEM自由度為28 800)的模型,僅耗時3.56 min,表明本文方法的計算效率明顯更高。另外,隨著頻率升高,文獻[4]的快速BEM會因為不支持高頻問題而失效,而本文的快速定向壓縮BEM仍能保證較高的計算效率。

觀察圖10和圖11還發(fā)現(xiàn),本文方法的計算結(jié)果在共振峰處與解析解存在較大誤差。對于圖10來說,這主要是因為其對應(yīng)的聲振耦合解析模型為一聲振耦合內(nèi)問題,結(jié)構(gòu)建模時沒有考慮阻尼,不存在能量耗散,即在共振峰處解析解理論上為無窮大,因此數(shù)值解和解析解在共振峰處存在一定誤差。在圖11中,主要是在一些比較尖銳的共振峰處數(shù)值解與解析解存在較大誤差,而在一些比較“鈍”的共振峰處(主要是在第1、2、3、5個共振峰處),二者吻合良好,這是因為該模型為一聲振耦合外問題,雖然結(jié)構(gòu)沒有阻尼,但在無窮遠處存在能量耗散,因此這些峰值處存在類似于阻尼耗能的現(xiàn)象[2]。當(dāng)理論計算結(jié)果趨于無窮時,數(shù)值方法的計算精度通常很難保證,這也是很多數(shù)值方法的局限性,如文獻[4,14]中的計算結(jié)果,也存在類似現(xiàn)象。

4.3 大規(guī)模工程應(yīng)用算例

在航海工程領(lǐng)域,聲對抗技術(shù)的發(fā)展對以潛艇、魚雷等為代表的水下航行器的振動和噪聲水平提出了更為嚴(yán)苛的標(biāo)準(zhǔn)。雖然水下結(jié)構(gòu)的剛度比航天器的大,但因聲介質(zhì)是水,其阻抗遠大于空氣阻抗,故聲振效應(yīng)尤為明顯。復(fù)雜殼體作為水下航行器的原型結(jié)構(gòu),對其聲振特性進行分析將能夠有效指導(dǎo)設(shè)計,降低研發(fā)成本。另外,此類問題屬于聲振耦合外問題,需要對外聲場進行建模,而這正是本文構(gòu)造的快速BEM-FEM耦合方法的優(yōu)勢所在。

建立如圖12所示殼體的半剖模型,軸向總長度l=6.0 m,半球端面半徑r=0.5 m,內(nèi)外殼壁厚t1=0.02 m,加強筋和隔板厚度t2=0.05 m。該系統(tǒng)的結(jié)構(gòu)和聲場的物理參數(shù)與前文模型1和模型2的相同。圖12中的箭頭表示幅值為1 N/m2的簡諧載荷,分別沿X方向和Z方向。該模型的結(jié)構(gòu)部分被離散成20 702個8節(jié)點四邊形單元,其中耦合界面上的單元數(shù)為6 900;外聲場被離散為13 800個6節(jié)點三角形單元;BEM-FEM系統(tǒng)的總自由度為420 654。

圖12 水中復(fù)雜殼體的半剖FEM模型

(a)X方向

(b)Y方向

(c)Z方向圖13 圖12中1、2和3位置處的位移頻響曲線

對該模型進行掃頻計算,關(guān)注的頻率范圍為1～1 000 Hz,步長取為1 Hz。提取圖12中1、2和3位置處的位移,分別繪制其X、Y和Z方向的頻響曲線,見圖13。提取圖13中第1個共振峰處(135 Hz)的解向量繪制位移云圖和聲壓云圖,見圖14。該算例共計耗時209 h,內(nèi)存占用13.8 GB。

(a)位移

(b)聲壓圖14 圖13頻響曲線第1個共振峰處(135 Hz)的響應(yīng)

5 總 結(jié)

本文采用邊界元-有限元耦合的方法求解結(jié)構(gòu)的聲振問題,并通過適用于寬頻波動問題的快速定向算法對常規(guī)邊界元進行加速,同時考慮耦合界面上聲場網(wǎng)格和結(jié)構(gòu)網(wǎng)格不匹配的情形,基于加權(quán)余量法建立了非匹配網(wǎng)格之間的數(shù)據(jù)交換矩陣,最終建立了一套性能優(yōu)異的聲振耦合求解方法。

(1)針對耦合界面上的非匹配網(wǎng)格數(shù)據(jù)交換問題,采用加權(quán)余量法構(gòu)造數(shù)據(jù)交換矩陣,并通過基于單元的積分方案計算耦合場上的單元積分。當(dāng)耦合界面為平面時,數(shù)據(jù)交換過程不損失精度,當(dāng)其為曲面時,數(shù)據(jù)交換的最大誤差不超過10-3。

(2)針對常規(guī)邊界元矩陣為滿陣的問題,采用適用于寬頻波動問題的快速定向算法對常規(guī)邊界元進行加速,實現(xiàn)了聲場的快速計算。與現(xiàn)有同類方法相比,本文方法的適用頻段范圍具有明顯優(yōu)勢,可達現(xiàn)有同類方法的3～4倍。

(3)將本文方法應(yīng)用于實際工程中常見的水下復(fù)雜殼體的聲振耦合外問題,在普通工作站上實現(xiàn)了自由度超過40萬的聲振耦合系統(tǒng)的響應(yīng)計算。