劉金實(shí)

(江蘇科技大學(xué) 船舶與海洋工程學(xué)院，江蘇 鎮(zhèn)江212003)

0 引 言

邊界元法是水下結(jié)構(gòu)振動(dòng)聲輻射預(yù)報(bào)的重要方法，與完全匹配層、無(wú)限單元法［1－3］相比，它具有直接求解開(kāi)放聲場(chǎng)、不需要在結(jié)構(gòu)外部添加額外網(wǎng)格的優(yōu)勢(shì)，然而由于其系統(tǒng)矩陣為密集矩陣，隨著問(wèn)題規(guī)模的增加，其計(jì)算耗時(shí)和內(nèi)存占用都迅速增長(zhǎng)，因此如何降低邊界元法的計(jì)算復(fù)雜度成了過(guò)去20年中學(xué)者研究的焦點(diǎn)，成果主要可分為兩大類:以提高矩陣－向量乘積運(yùn)算速度為目標(biāo)的算法和矩陣壓縮算法。前者以快速多極子邊界元法［4］為代表，后者更接近于純代數(shù)方法，以層級(jí)矩陣技術(shù)和自適應(yīng)交叉近似算法為代表。由于快速多極子算法的程序?qū)崿F(xiàn)受單元階次的影響較大，多極子樹(shù)的構(gòu)建依賴于模型的幾何特征，在工程應(yīng)用中難以推廣。

層級(jí)矩陣的核心思想是對(duì)密集矩陣進(jìn)行層層分割，將子塊的位置、大小等信息按樹(shù)狀的數(shù)據(jù)結(jié)構(gòu)加以保存，對(duì)于其中秩較低的子塊進(jìn)行奇異值分解(SVD)，保留其余分塊的密集格式。在之后的數(shù)年中，該方法被應(yīng)用于聲場(chǎng)、電磁場(chǎng)問(wèn)題的計(jì)算中。在層級(jí)矩陣的基礎(chǔ)上，Bebendorf［5］于2004年提出了自適應(yīng)交叉近似(ACA)算法，該算法主要從2個(gè)方面完善并改進(jìn)了層級(jí)矩陣法，首先是利用交叉近似方法對(duì)矩陣中的低秩子塊進(jìn)行逼近，替代了SVD算法，從而降低了對(duì)矩陣進(jìn)行壓縮的計(jì)算復(fù)雜度，其次是通過(guò)考察積分核的退化特性使算法具備了自適應(yīng)的能力。Grasedyck［6］在交叉逼近的過(guò)程中引入了參考行和參考列以改進(jìn)ACA 算法的逼近策略，避免了原有算法在少數(shù)情況下崩潰的風(fēng)險(xiǎn)，也使得ACA 算法的收斂速度有了輕微提高，同時(shí)首次將ACA 算法應(yīng)用于邊界元的快速計(jì)算，形成了ACA－BEM 方法。鄭昌軍等［7］將ACA－BEM 應(yīng)用于聲靈敏度的計(jì)算中，實(shí)現(xiàn)了聲場(chǎng)問(wèn)題最優(yōu)化的快速計(jì)算。

1 理論基礎(chǔ)

1.1 低秩近似理論

假設(shè)一個(gè)秩為k的矩陣M∈Cm×n，M的R(k)矩陣低秩近似可用因式分解的形式表示為:

其中U∈Cm×k;V∈Cn×k。式(1)也可更方便地表示為:

經(jīng)過(guò)式(1)的分解，存儲(chǔ)M 所需的內(nèi)存空間由O(m × n)變?yōu)镺(k × (m + n))，設(shè)有任意向量t∈Cn×1，則利用式(2)計(jì)算M × t，所需的浮點(diǎn)數(shù)乘法操作次數(shù)也由O(m × n)變?yōu)镺(k × (m + n))。如圖1所示，當(dāng)M的秩遠(yuǎn)小于其維數(shù)時(shí)，利用式(1)可以顯著降低其存儲(chǔ)空間，并提高M(jìn) × t的運(yùn)算效率。

圖1 矩陣分解示意圖Fig.1 Diagrammatic sketch of matrix decomposition

然而在實(shí)際應(yīng)用中，精確求出M的秩需要較大的計(jì)算量，會(huì)使得近似失去意義，為此須引入ε－秩及最佳近似概念:

其中M的最佳近似矩陣N∈Cm×n;rank(N)=k;‖…‖代表矩陣的范數(shù)。

1.2 自適應(yīng)交叉近似算法

現(xiàn)有非零矩陣M∈Cm×n，且其所有元素的值均為已知。全選主元交叉近似算法首先利用以下步驟找出矩陣M的秩為1的近似矩陣M1:

1)確定M 中模值最大的元素所在行、列，記為i1和j1，并設(shè)δ=Miljl;

2)取出矩陣M 中第i1行和第j1列的向量，存儲(chǔ)矩陣U*的第一列u1=M，j1，矩陣V*的第一行v1=Mi1，/δ。

則矩陣M1=U*V*。設(shè)整數(shù)p:1 ≤p ≤k，設(shè)第p步逼近完成后的余項(xiàng)為Rp，則將Rp作為第p+1步逼近的目標(biāo)，從而形成如式(4)所示的關(guān)系。

在每一步逼近結(jié)束時(shí)，都可利用式(5)考察誤差閾值是否已得到滿足。

1.3 自適應(yīng)交叉近似算法

準(zhǔn)備階段:輸入矩陣M∈Cm×n，記‖M‖2= τ，設(shè)定收斂誤差閾值ε，定義近似矩陣的Frobenius 范數(shù)為F0，隨機(jī)生成整數(shù)il，滿足1≤il≤m，設(shè)已使用的行坐標(biāo)集合D={i1};

1)設(shè)定迭代計(jì)數(shù)p = 1，計(jì)算矩陣M的第i1行向量Mi1，*，確定Mi1，*中模值最大的列坐標(biāo)j1，即j1=argmax(Mi1，*)，以及該元素的值δ=Mi1，j1;

2)計(jì)算出M的第j1列，記為M*，j1，根據(jù)式(6)得到逼近向量u1和v1:

3)利用式(7)更新逼近矩陣的Frobenius 范數(shù):

4)利用式(8)判別逼近是否已經(jīng)滿足誤差閾值的要求，若滿足則跳轉(zhuǎn)至步驟8)，否則繼續(xù)執(zhí)行步驟5)。

5)D=D ∪{ip}，根據(jù)列向量up從{1，…，m}/D 中選取行坐標(biāo)ip+1，即ip+1=argmax(up)，計(jì)算Mip+1，*，利用式(9)計(jì)算余項(xiàng)矩陣第ip+1行作為逼近向量vp+1。

6)令jp+1=argmax(vi+1)，根據(jù)式(10)計(jì)算R*，jp+1，取出交叉點(diǎn)的余項(xiàng)值δ=Rip+1，jp+1，若δ 則終止循環(huán)，跳轉(zhuǎn)至8)，否則利用式(11)計(jì)算逼近向量up+1。

7)利用式(12)更新逼近矩陣的Frobenius 范數(shù)并跳轉(zhuǎn)至步驟4)。

8)利用式(13)給出矩陣U*∈Cm×p和V*∈Cp×n，算法終止。

圖2 交叉近似示意圖Fig.2 Diagrammatic sketch of cross approximation

1.4 基于主成分分析(PCA)的求解域分割

由于邊界元矩陣中的行、列坐標(biāo)都對(duì)應(yīng)了網(wǎng)格中的各個(gè)節(jié)點(diǎn)，而矩陣的子塊則代表著2個(gè)節(jié)點(diǎn)子集之間的相互作用，因此要將ACA 算法應(yīng)用于對(duì)邊界元矩陣的壓縮，首先應(yīng)根據(jù)幾何相容性條件對(duì)求解域進(jìn)行分割，得到若干節(jié)點(diǎn)集，通過(guò)將節(jié)點(diǎn)集之間的相互作用映射為邊界元矩陣中的分塊，利用幾何相容性條件即可對(duì)子塊的低秩特性進(jìn)行預(yù)判。

文獻(xiàn)［8－9］中較為常用的笛卡爾坐標(biāo)系分割法，其思路較為簡(jiǎn)單，即首先找到包圍求解域的最小長(zhǎng)方體，長(zhǎng)方體的各個(gè)面均與笛卡爾坐標(biāo)系的坐標(biāo)軸垂直，然后沿著3個(gè)坐標(biāo)軸對(duì)長(zhǎng)方體進(jìn)行多級(jí)等分，直到所得到的子區(qū)域足夠小或者滿足幾何相容條件。當(dāng)求解域?yàn)槎S空間中的一個(gè)圓時(shí)，分割結(jié)果如圖3所示。

圖3 圓形求解域分割示意圖Fig.3 Schematic diagram of splitting circular solving domain

根據(jù)積分核退化理論［5］，設(shè)幾何相容參數(shù)η=，則由圖3 可看出，當(dāng)包圍求解域的長(zhǎng)方體接近于正方體時(shí)，利用笛卡爾坐標(biāo)系分割法得到的子塊之間幾乎都滿足幾何相容條件。然而當(dāng)求解域?yàn)橐粋€(gè)長(zhǎng)寬比為4的長(zhǎng)方形時(shí)，笛卡爾坐標(biāo)系分割結(jié)果如圖4(a)所示，其中a與c、b與d 之間均不滿足幾何相容條件，但如果只沿長(zhǎng)方形的長(zhǎng)邊進(jìn)行分割，如圖4(b)所示，則任意2個(gè)不同子塊之間都滿足幾何相容條件。由此可見(jiàn)，對(duì)求解域的分割方式影響著ACA 算法的壓縮效果，進(jìn)而可能影響到ACA 邊界元法整體的計(jì)算效率。

實(shí)際分析中考察的模型既可能是如圖3和圖4 一樣的形狀規(guī)則，更多情況下也可能是各種不規(guī)則的形狀，因此必須通過(guò)找到一種適應(yīng)不同形狀的求解域分割算法來(lái)提高子塊間滿足幾何相容條件的概率。主 成 分 分 析［10］(Principle Component Analysis，PCA)又被稱為主分量分析。利用主成分分析對(duì)數(shù)據(jù)進(jìn)行處理，可以從影響問(wèn)題的多個(gè)變量中選出影響較大的少數(shù)變量，這種思想在圖像分析甚至經(jīng)濟(jì)學(xué)等領(lǐng)域有著廣泛的應(yīng)用。為了說(shuō)明PCA 在分割邊界元求解域中的應(yīng)用，首先考慮一個(gè)二維空間中點(diǎn)集z，包含的節(jié)點(diǎn)數(shù)為n，節(jié)點(diǎn)分布如圖5所示，其中節(jié)點(diǎn)i∈［1，n］的坐標(biāo)記為zi，z的幾何中心坐標(biāo)為Cz，則z的主方向w的定義由下式給出:

圖4 橫向分割示意圖和笛卡爾坐標(biāo)系分割示意圖Fig.4 Subdivision of horizontal splitting and subdivision of cartesian splitting

其中v為任意單位向量。在實(shí)際計(jì)算中可通過(guò)計(jì)算ZTZ最大特征值所對(duì)應(yīng)的特征向量得到主方向向量w，即:

圖5 笛卡爾坐標(biāo)系分割和主成分分析法分割Fig.5 Subdivision of cartesian splitting and subdivision of PCA splitting

將w 作為坐標(biāo)軸，得到z 中各點(diǎn)的坐標(biāo)，在坐標(biāo)值的最大、最小值的中心處取與w 垂直的分割面d，將z 分割為兩部分，如圖5(b)所示。

在對(duì)網(wǎng)格進(jìn)行逐級(jí)分割的過(guò)程中，需要重復(fù)計(jì)算ZTZ 及其特征值、特征向量，相對(duì)于傳統(tǒng)方法引入了額外的計(jì)算量，但由于在同一模型的各個(gè)分析工況間并不需要重復(fù)運(yùn)算，因此不會(huì)對(duì)邊界元分析的整體計(jì)算復(fù)雜度造成影響。

2 數(shù)值仿真

2.1 算法驗(yàn)證

為了驗(yàn)證ACA 邊界元法的有效性，分析算法各參數(shù)設(shè)置對(duì)計(jì)算精度與效率的影響以及迭代求解器的收斂性能，本節(jié)采用圓柱殼作為驗(yàn)證模型，其中球殼半徑為1 m，圓柱殼長(zhǎng)為2 m，半徑0.3 m，模型外部介質(zhì)聲速為1 500 m/s，密度為1 000 kg/m3。模型的形心坐標(biāo)為(0，0，0)，且圓柱殼的軸向與z軸平行。統(tǒng)一設(shè)定求解域分割的最小子塊含節(jié)點(diǎn)數(shù)不低于100，誤差參數(shù)ε=10－6。

利用映射網(wǎng)格劃分方法，將模型劃分為746 節(jié)點(diǎn)、744個(gè)四邊形單元的網(wǎng)格，設(shè)定所有節(jié)點(diǎn)沿外法線方向的振速為1，在(10，0，0)出設(shè)置一個(gè)聲壓觀測(cè)點(diǎn)。分別利用傳統(tǒng)邊界元法和ACA 邊界元法計(jì)算聲壓觀測(cè)點(diǎn)處的輻射聲壓級(jí)，結(jié)果的誤差曲線如圖6所示。

圖6 圓柱殼模型輻射聲壓級(jí)對(duì)比圖Fig.6 Comparison of radiate sound pressure level for cylindrical shell

由圖6 可看出，當(dāng)誤差控制參數(shù)設(shè)定為ε=10－6時(shí)，ACA－BEM與傳統(tǒng)邊界元法的預(yù)報(bào)結(jié)果具有很好的一致性。

2.2 計(jì)算效率分析

除了2.1 節(jié)中的圓柱殼模型，本節(jié)還考慮一個(gè)半徑為1 m、厚度為0.003 m的水下球殼模型，將2個(gè)模型分別劃分8 組不同密度的網(wǎng)格，如表1所示。令模型表面所有節(jié)點(diǎn)的振速為1 m/s，求解器內(nèi)存占用和計(jì)算耗時(shí)曲線如圖7和圖8所示。

表1 網(wǎng)格規(guī)模表Tab.1 Table of mesh sizes

圖7 內(nèi)存占用量對(duì)比圖Fig.7 Comparison of BEM matrices memory cost

圖8 求解耗時(shí)對(duì)比圖Fig.8 Comparison of solving time cost

由圖7和圖8 可看出，由于采用了迭代求解器，ACA 邊界元法的求解器只需要額外存儲(chǔ)預(yù)條件矩陣的稀疏LU 分解因子，而對(duì)密集矩陣的LU 分解，分解因子和原矩陣本身的規(guī)模基本一致，因此在求解器的內(nèi)存占用量上，ACA 邊界元法具有較大的優(yōu)勢(shì);從增長(zhǎng)趨勢(shì)來(lái)看，ACA 邊界元法的求解運(yùn)算耗時(shí)隨節(jié)點(diǎn)數(shù)的增長(zhǎng)速度顯著緩于傳統(tǒng)邊界元法，但在網(wǎng)格節(jié)點(diǎn)數(shù)小于2 000 時(shí)并未顯現(xiàn)出優(yōu)勢(shì)，主要原因是傳統(tǒng)邊界元法的求解過(guò)程中計(jì)算量最大的LU分解采用了LAPACK 標(biāo)準(zhǔn)函數(shù)，運(yùn)行時(shí)的調(diào)度優(yōu)化水平很高，而迭代求解器中矩陣與向量的乘法操作被替代為大量小型矩陣與向量相乘，雖然每次小矩陣與向量相乘都采用了優(yōu)化程度較高的BLAS 標(biāo)準(zhǔn)子程序，能夠調(diào)動(dòng)CPU 所有計(jì)算核心并行運(yùn)算，但每個(gè)線程的運(yùn)行時(shí)間都非常短，頻繁的線程切換浪費(fèi)了大量的時(shí)間，據(jù)筆者觀察，在計(jì)算程序運(yùn)行過(guò)程中，CPU的峰值占用量不足35%。因此，通過(guò)改進(jìn)計(jì)算程序的并行方式，充分利用機(jī)器資源，還可以較大幅度的提高計(jì)算程序的效率。

3 結(jié) 語(yǔ)

1)采用ACA－BEM 算法替代傳統(tǒng)邊界元法，并利用主成分分析法分割求解域，建立圓柱殼體模型，利用均勻脈動(dòng)殼體問(wèn)題驗(yàn)證了方法的有效性和準(zhǔn)確性。

2)對(duì)比了ACA－BEM和傳統(tǒng)邊界元法在不同模型尺寸比、不同網(wǎng)格規(guī)模條件下的計(jì)算效率和內(nèi)存占用量。ACA－BEM的計(jì)算效率受結(jié)構(gòu)長(zhǎng)寬比的影響較小，當(dāng)節(jié)點(diǎn)數(shù)不足2000 時(shí)，ACA－BEM的計(jì)算效率與傳統(tǒng)邊界元相比不具備優(yōu)勢(shì)，但隨著節(jié)點(diǎn)數(shù)的繼續(xù)增多，ACA－BEM 在計(jì)算效率和內(nèi)存占用量?jī)煞矫娑硷@現(xiàn)出明顯的優(yōu)勢(shì)。

［1］JEAN P B.A perfectly matched layer for the absorption of electromagnetic waves ［J］.Journal of Computational Physics，1994，114(2):185－200.

［2］LUC C，F(xiàn)YFE K R，COYETTE J P.A variable order infinite acoustic wave envelope element［J］.Journal of Sound and Vibration，1994，171(4):483－508.

［3］BURNETT D S.A three-dimensional acoustic infinite element based on a prolate spheroidal multipole expansion［J］.Journal of Acoustic Society of America，1996，96:2798－2816.

［4］MICHAEL－A E，DEMBART B.Multipole translation theory for the three-dimensional laplace and helmholtz equations［J］.SIAM Journal on Scientific Computing，1995，16(4):865－897.

［5］MARIO B，RJASANOW S.Adaptive low-rank approximation of collocation matrices［J］.Computing，2003，70(1):1－24.

［6］LARS G.Adaptive recompression of－matrices for BEM［J］.Computing，2005，74(3):205－223.

［7］鄭昌軍，陳磊磊，陳海波.基于自適應(yīng)交叉近似邊界元法的快速聲學(xué)靈敏度分析［J］.計(jì)算力學(xué)學(xué)報(bào)，2014(4):413－418.

［8］吳君輝，曹祥玉，袁浩波，等.自適應(yīng)交叉近似算法在矩量法中的應(yīng)用［J］.空軍工程大學(xué)學(xué)報(bào)(自然科學(xué)版)，2013(5):76－79.

［9］楊劍.半空間環(huán)境中的自適應(yīng)交叉近似方法研究［D］.西安:西安電子科技大學(xué)，2013.

［10］Ian Jolliffe.Principal component analysis［M］.Wiley Online Library，2005.