不連續(xù)伽遼金積分方程區(qū)域分解方法在定常迭代中的收斂性研究

(1.成都師范學(xué)院物理與工程技術(shù)學(xué)院 四川 成都 611130；2.電子科技大學(xué)能源科學(xué)與工程學(xué)院 四川 成都 611731)

引言

表面積分方程方法在電磁散射和輻射問題中得到了廣泛的應(yīng)用。為了降低矩量法求解積分方程時所面臨的O(N2)的計算復(fù)雜度和存儲復(fù)雜度，快速算法獲得了極大地發(fā)展[1]。但對于電大問題來說，其求解過程仍需要大量的計算機(jī)存儲資源，而這可能超出當(dāng)前平臺的計算能力。為了降低電大問題求解過程中對計算機(jī)內(nèi)存資源的需求，提高當(dāng)前平臺的計算能力，可以通過逐個求解子問題的方式來求解原問題。為此，采用經(jīng)典的定常迭代思路來求解原始問題的區(qū)域分解方法應(yīng)運(yùn)而生。文獻(xiàn)[2]提出了一種前后向迭代方案，引入了一定長度的緩沖區(qū)以抑制不必要的偽邊緣效應(yīng)。文獻(xiàn)[3]采用了多邊緩沖區(qū)，實現(xiàn)了更快的收斂速度。文獻(xiàn)[4]使用了一層三角形單元作為緩沖區(qū)，降低了子域問題中緩沖區(qū)上的額外未知量數(shù)目，并通過引入電流連續(xù)性條件的方式，抑制偽邊緣效應(yīng)，獲得了穩(wěn)定的收斂效果。

以上方法均屬于重疊型區(qū)域分解方法。文獻(xiàn)[5]提出了一種非重疊型區(qū)域分解方法，可以采用定常迭代法進(jìn)行求解，以實現(xiàn)較少的存儲資源占用。為了在不引入人工端面的條件下，構(gòu)造非重疊型區(qū)域分解模型，并使用定常迭代法進(jìn)行求解，文獻(xiàn)采用了顯示邊界條件來約束邊界處兩個相鄰的半Rao-Wilton-Glisson(RWG)基函數(shù)的系數(shù)，獲得與原問題等價的區(qū)域分解問題，同時生成一個超定的線性方程系統(tǒng)。該超定線性系統(tǒng)通常需要法方程方法進(jìn)行求解，這可能導(dǎo)致迭代收斂相對變慢。最近，文獻(xiàn)提出了一種基于不連續(xù)伽遼金方法的非重疊型區(qū)域分解方法，并采用Krylov子空間迭代法進(jìn)行求解。但當(dāng)采用定常迭代法進(jìn)行求解時，由于其引入的內(nèi)罰穩(wěn)定參數(shù)具有明顯的頻率相關(guān)性，導(dǎo)致外層迭代的收斂性對工作頻率十分敏感，甚至在某些頻段內(nèi)收斂緩慢。

為了在寬頻帶內(nèi)實現(xiàn)定常迭代法的快速收斂，本文使用了一個與頻率無關(guān)的內(nèi)罰穩(wěn)定參數(shù)，用于求解采用均勻網(wǎng)格剖分的電大導(dǎo)體目標(biāo)的電磁散射問題。這個參數(shù)正比于平均剖分尺寸，并反比于工作波長。為了驗證該參數(shù)對迭代收斂性的影響，對給定了幾何尺寸的某一目標(biāo)，本文首先考察了該頻率無關(guān)的穩(wěn)定參數(shù)和原始參數(shù)，在頻率范圍為0.1GHz到10GHz下的迭代矩陣譜半徑。然后對比考察了相同工作頻率下，使用這兩個參數(shù)的迭代收斂速度。此外，本文還討論了不同子域數(shù)目對收斂性和精度的影響。

一、區(qū)域分解方法回顧

考慮理想導(dǎo)體目標(biāo)在自由空間中的散射問題。假設(shè)目標(biāo)表面被分為M個子域S=S1∪…∪SM。子域內(nèi)部的感應(yīng)電流采用RWG基函數(shù)展開，跨邊界電流則采用單極RWG基函數(shù)展開。利用導(dǎo)體表面的切向邊界條件，并采用伽遼金測試方法，可以建立如下的電場積分方程(EFIE)和磁場積分方程(MFIE)：

(1)

(2)

其中

(3)

(4)

(5)

CFIE=αEFIE+(1-α)η0×MFIE,α∈(0,1)

(6)

式中，α為比例因子，通常取0.5。

二、定常迭代形式

對上述方程采用分域基函數(shù)離散，并采用伽遼金方法測試后，可以獲得如下的線性系統(tǒng)(假設(shè)M=2)：

(7)

其中Am是子域Sm上的CFIE阻抗矩陣，Bm,n是子域Sm到Sn的耦合矩陣。為了利用有限的內(nèi)存資源求解電大問題，上述矩陣方程可以采用定常迭代法來求解，允許計算機(jī)一次僅處理一個子域問題。該迭代過程可以表示為：

(8)

對于式(8)所對應(yīng)的子域Sm上的矩陣方程，可以采用 Krylov子空間迭代法GMRES進(jìn)行求解。求解式(8)的一次迭代稱作內(nèi)層迭代，而求解所有子域的一次迭代稱為外層迭代。在第k外層迭代后的相對殘差定義為：

(9)

其中‖·‖表示一個復(fù)向量的二范數(shù)。當(dāng)ε(k)小于給定的誤差門限后，外層迭代停止。多層快速多極子方法[3]被用來加速(8)中的矩矢相乘運(yùn)算。

三、收斂性分析

Mx(k+1)=Nx(k)+b

(10)

其中A=M-N是對原始矩陣A的一個分裂。在本文中，分裂后的矩陣M和N可以具體表示為：

(11)

(10)式是否收斂到真解x=A-1b，取決于迭代矩陣G=M-1N的特征值。如果G的譜半徑ρ(G)=max{|λ|:λ∈λ(G)}小于1.0，則由(10)式定義的近似解序列{x(k)}將會收斂于x=A-1b。此外，譜半徑越小，收斂速度越快。

圖1 金屬球的表面分區(qū)

然后我們考察了固定目標(biāo)的幾何尺寸，改變工作頻率對定常迭代收斂性的影響。在這個算例中，目標(biāo)的電尺寸是隨著頻率的變化而變化的。考察一個半徑為r=1m的金屬球，工作頻率分別選取為0.1GHz和3GHz。同樣地，該目標(biāo)被分為兩個子域，如圖1所示，并采用邊長為0.1λ0的三角形網(wǎng)格進(jìn)行均勻剖分。表1給出了采用不同穩(wěn)定參數(shù)收斂到ε=10-3時的外層迭代次數(shù)?？梢钥吹剑捎帽疚膮?shù)時，定常迭代的收斂性不隨頻率的變化而變化，且保持較快的收斂速度。

圖2 不同穩(wěn)定參數(shù)下的迭代矩陣譜半徑

表1不同穩(wěn)定參數(shù)下的迭代收斂情況(ε=10-3)

Freq.(GHz)0.11β=0.15|loghλ-10|77β=0.1|logh|157

四、數(shù)值算例

在這一部分，數(shù)值算例用于驗證改進(jìn)方法的正確性和有效性。本文數(shù)值算例在一個配置為雙核CPU(頻率為3.2GHz)和2GB內(nèi)存，搭載Win64操作系統(tǒng)的個人計算機(jī)上完成。算例1為金屬立方體目標(biāo)，考察不同子域數(shù)目對外層迭代收斂性的影響。算例2為金屬球為目標(biāo)，用于驗證分區(qū)數(shù)目對計算精度的影響。

圖3 不同子域數(shù)目下的收斂情況

圖4 目標(biāo)區(qū)域分解示意圖

圖5 目標(biāo)計算結(jié)果對比

算例1 選擇2m×2m×2m的金屬立方體作為計算目標(biāo)。將該立方體分別劃分為 4、16和56個子域，如圖3所示。入射波頻率為300MHz，入射角度為 θ=0°，φ=0°，電場沿θ方向極化。圖3 給出了不同分區(qū)數(shù)目下該方法的迭代收斂曲線。從4個子域到 56 個子域，收斂到 10-2所需的迭代次數(shù)從 8 增加到 11；而收斂到 10-6所需的迭代次數(shù)從 26 增加到 34。從該算例可以看到，在子域數(shù)目增加的情況下，該方法穩(wěn)健性良好。此外，采用4,16,56個分區(qū)所需計算內(nèi)存和所需計算時間跟不采用區(qū)域分解方法所需計算內(nèi)存和計算時間對比可知，使用區(qū)域分解方法能有效的降低求解過程中的峰值內(nèi)存需求，但與此同時，計算時間會有所增加。

算例2 考慮半徑為 1m 的金屬球，如圖4 所示，在 300MHz 平面波照射下的電磁散射。目標(biāo)分別8 個子域。為了獲得更多的分區(qū)，采用邊長為 0.1 波長的三角形單元對原目標(biāo)進(jìn)行剖分，并將每個單元視為一個區(qū)域，共劃分出 2586 個子域，平面波入射角度為 θ=0°，φ=0°，電場沿θ方向極化。對于一般的幾何模型，改進(jìn)后的方法在提高單機(jī)求解能力的基礎(chǔ)上，較原方法能進(jìn)一步降低求解時間，提高了求解效率。圖5給出了使用兩種穩(wěn)定參數(shù)下的計算結(jié)果對比，吻合良好。

五、總結(jié)

使用與頻率無關(guān)的穩(wěn)定參數(shù)，并選取合適的系數(shù)，使定常迭代法在寬頻帶內(nèi)穩(wěn)定、快速收斂，可在降低內(nèi)存需求的基礎(chǔ)上，進(jìn)一步降低求解時間，提高當(dāng)前計算平臺的求解能力。在固定目標(biāo)電尺寸和幾何尺寸的情況下，使用本文給出的穩(wěn)定參數(shù)，定常迭代的收斂性均不再表現(xiàn)出與工作頻率的相關(guān)性。