劉金實

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

0 引 言

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

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

1 理論基礎(chǔ)

1.1 低秩近似理論

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

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

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

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

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

其中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列的向量，存儲矩陣U*的第一列u1=M，j1，矩陣V*的第一行v1=Mi1，/δ。

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

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

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

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

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

2)計算出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 中選取行坐標ip+1，即ip+1=argmax(up)，計算Mip+1，*，利用式(9)計算余項矩陣第ip+1行作為逼近向量vp+1。

6)令jp+1=argmax(vi+1)，根據(jù)式(10)計算R*，jp+1，取出交叉點的余項值δ=Rip+1，jp+1，若δ 則終止循環(huán)，跳轉(zhuǎn)至8)，否則利用式(11)計算逼近向量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)的求解域分割

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

文獻［8－9］中較為常用的笛卡爾坐標系分割法，其思路較為簡單，即首先找到包圍求解域的最小長方體，長方體的各個面均與笛卡爾坐標系的坐標軸垂直，然后沿著3個坐標軸對長方體進行多級等分，直到所得到的子區(qū)域足夠小或者滿足幾何相容條件。當求解域為二維空間中的一個圓時，分割結(jié)果如圖3所示。

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

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

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

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

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

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

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

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

2 數(shù)值仿真

2.1 算法驗證

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

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

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

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

2.2 計算效率分析

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

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

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

圖8 求解耗時對比圖Fig.8 Comparison of solving time cost

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

3 結(jié) 語

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

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

［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)交叉近似邊界元法的快速聲學靈敏度分析［J］.計算力學學報，2014(4):413－418.

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

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

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