刁愛民，楊慶超，王杏青

(海軍工程大學(xué) a.科研部;b.電氣與信息工程學(xué)院 電磁環(huán)境與防護(hù)系，武漢 430033)

海水中電偶極子的電磁場(chǎng)是目前水下通訊和水下目標(biāo)探測(cè)方面的研究重點(diǎn)。電偶極子是指由離開一定距離的電量相同、符號(hào)相反的兩個(gè)點(diǎn)電荷組成的系統(tǒng)，兩個(gè)點(diǎn)電荷之間的距離與待求場(chǎng)點(diǎn)到他們中心的距離相比很小，是最簡(jiǎn)單也是最基本的一種電磁場(chǎng)場(chǎng)源，求解電偶極子產(chǎn)生的電場(chǎng)也就可計(jì)算出任意場(chǎng)點(diǎn)的電場(chǎng)值。

海水是一種導(dǎo)電媒質(zhì)，當(dāng)海水的深度較大時(shí)，假定海水的電導(dǎo)率恒定不變，可將海水看作無限大導(dǎo)體空間，并利用有限元方法對(duì)感興趣的區(qū)域進(jìn)行有限分割，選取適當(dāng)?shù)牟逯岛瘮?shù)，就可近似推導(dǎo)電場(chǎng)的解析表達(dá)式。當(dāng)海水深度較小時(shí)，可將深海等效為空氣-海水兩層模型。

1 計(jì)算模型

1.1 深海點(diǎn)電荷模型［1］

將海水的電導(dǎo)率視為恒定不變，如圖1 所示，電荷量為Q 的點(diǎn)電荷的電場(chǎng)公式為:

用一個(gè)半徑為r 的實(shí)心等勢(shì)鐵球體代替點(diǎn)電荷，若在球體上加上均勻的固定電勢(shì)φ0，則由高斯定律可得:

取另一半徑為R 的高斯面，由高斯定理得:

由E=εD，可得即電勢(shì)為φ0，半徑為r 的等勢(shì)球體與帶電量為Q=4πεrφ0的點(diǎn)電荷產(chǎn)生的遠(yuǎn)場(chǎng)一樣，可以相互代替。即深海點(diǎn)電荷模型可以是一個(gè)等勢(shì)小球。

圖1 點(diǎn)電荷模型

2 深海中電偶極子有限元模型

由點(diǎn)電荷模型可知，可用2 個(gè)實(shí)心等勢(shì)鐵球體來代替點(diǎn)電荷構(gòu)建一個(gè)電偶極子模型，如圖2 所示。由于電偶極子是指由相距很小距離的電量相同，符號(hào)相反的2 個(gè)點(diǎn)電荷組成的系統(tǒng)，故應(yīng)滿足2 個(gè)條件:①小球的半徑比二者之間的距離至少小1 個(gè)數(shù)量級(jí);②求解的范圍比二者之間的距離至少高1 個(gè)數(shù)量級(jí)。即該模型中電場(chǎng)源為2 個(gè)電勢(shì)相反的球體，求解域?yàn)槌浞执蟮某錆M海水的立方體，邊界條件為自然邊界條件。

圖2 電偶極子有限元模型

2 電偶極子電場(chǎng)的微分方程和邊界條件［2］

2.1 深海電偶極子有限元模型

式中:φ 表示介質(zhì)中以及金屬表面的電位值;φ0是邊界上電位的給定值，該模型中邊界的電位為0。

2.2 空氣-海水兩層有限元模型［3］

當(dāng)空間分成如圖3 所示的2 個(gè)平面，其中一半充滿空氣，另一半為不流動(dòng)的海水，空氣為區(qū)域1，海水為區(qū)域2。以空氣-海水界面為z =0 的平面，空氣占據(jù)z ＞0 的上半空間，而海水則占據(jù)下半空間，空氣的電磁特性為常數(shù)ε0、μ0和σ0(σ0=0)，假定海水的電磁特性為常數(shù)ε、μ0和σ，則兩部分介質(zhì)中的傳播常數(shù)分別為:

圖3 空氣—海水兩層模型

假設(shè)在海水中有一電流源J，取時(shí)諧因子為ejωt，則電場(chǎng)波動(dòng)方及邊界條件為［4］:

選取一個(gè)足夠大的六面體區(qū)域包圍電流源J，則在六面體的外邊界上，電磁場(chǎng)的邊界條件是:

3 有限元方程的求解

應(yīng)用廣義變分原理，對(duì)于研究的電磁場(chǎng)方程及對(duì)應(yīng)的邊界條件的變分為［5］:

采用的基本單元是六面體單元，并采用了與常規(guī)有限元不同的網(wǎng)格剖分技術(shù)。具體方法是將場(chǎng)域劃分成若干個(gè)網(wǎng)曲面，每個(gè)曲面層由若干個(gè)六面體單元組成。由于在海水平面剖分單元的密度要求高，而在所取區(qū)域上部空氣中和深部單元的密度要求低，因而在實(shí)際剖分時(shí)對(duì)于空氣中和不同深度海水中的曲面層是永不等距的網(wǎng)格剖分，分別采用向上收縮網(wǎng)格和向下收縮網(wǎng)格的方法，大大地減少了單元和節(jié)點(diǎn)數(shù)。經(jīng)求導(dǎo)，并整理，綜合后得總剛度矩陣K:

式中:K 為n×n 方陣;b 為1 ×n 列陣;n 為剖分節(jié)點(diǎn)數(shù)。

將其代入式(8)即可求得垂直電偶極子的電場(chǎng)有限元方程組。

將其代入式(8)即可求得水平電偶極子的電場(chǎng)有限元方程組。

對(duì)于式(8)的有限元方程，可通過多種解方程組的方法進(jìn)行求解，直接解法中有LU 分解、Gasuss 法和帶狀矩陣算法等;迭代法有共軛梯度法和復(fù)雙共軛梯度法。當(dāng)線性方程組系數(shù)矩陣的條件數(shù)較差時(shí)，LU 分解方法最為穩(wěn)定。如果系數(shù)矩陣中將散度效正項(xiàng)考慮進(jìn)去，則方程組的系數(shù)矩陣變?yōu)椴粚?duì)稱，系數(shù)矩陣需要存儲(chǔ)滿陣，這樣就很難滿足三維電磁問題的求解。復(fù)雙共軛梯度法收斂速度快，適用于大型稀疏對(duì)稱矩陣，采用一維壓縮內(nèi)存存儲(chǔ)方法，在計(jì)算機(jī)上可實(shí)現(xiàn)三維電磁問題的求解。

4 誤差分析

在深海模型中，海水的介電常數(shù)為81，球體半徑為0.1 mm，兩球間距離l=10 mm，球表面的電位為φ0= ±105V，則該模型等效于一個(gè)電荷量q=320π(C)，距離l=10 mm 的電偶極子產(chǎn)生的場(chǎng)［6］。

用Matlab 計(jì)算得到該電偶極子的場(chǎng)如圖4 所示。

圖4 電偶極子電勢(shì)分布

穿過電偶極子的對(duì)角線上的電場(chǎng)分布如圖5 所示。

圖5 電偶極子對(duì)角線上的電場(chǎng)分布

有限元法求得的兩小球的產(chǎn)生的電場(chǎng)分布如圖6所示［7］。

圖6 ansoft 仿真電偶極子的電場(chǎng)分布

將二者對(duì)角線的(-3，3)部分的電場(chǎng)分布進(jìn)行比較，結(jié)果如圖7 所示，有限元的計(jì)算結(jié)果變化較為平緩，誤差較大，原因是小球的半徑再小也不能近似為零，所以近場(chǎng)的誤差較大，符合實(shí)際。結(jié)果表明小球的曲率越大，附近的電場(chǎng)也就越大，計(jì)算結(jié)果也會(huì)更加陡峭，跟理論的結(jié)果更相近，但對(duì)剖分精度的要求加大，同時(shí)增加了計(jì)算量。

圖7 對(duì)角線上(-3，3)部分的電場(chǎng)分布比較

將二者對(duì)角線(-60，-8)部分的電場(chǎng)分布進(jìn)行比較，如圖8 所示，結(jié)果表明二者幾乎重合，表明小球的電偶極子模型的遠(yuǎn)場(chǎng)效果很好，同時(shí)對(duì)小球的半徑的數(shù)量級(jí)要求不高，較為容易實(shí)現(xiàn)。

圖8 對(duì)角線上(-60，-8)部分的電場(chǎng)分布比較

5 結(jié)束語(yǔ)

本文利用高斯定律從理論上分析了等勢(shì)小球和點(diǎn)電荷的電場(chǎng)分布的一致性，同時(shí)利用有限元方法推導(dǎo)出了海水中靜態(tài)和時(shí)諧電偶極子的電場(chǎng)有限元方程，通過對(duì)有限元方程的求解計(jì)算得到海水中任意場(chǎng)點(diǎn)的電場(chǎng)值。同時(shí)利用Matlab 建立的電偶極子模型比較、分析了靜態(tài)電偶極子模型的誤差，結(jié)果表明有限元電偶極子模型具有很好的遠(yuǎn)場(chǎng)效果。可見，運(yùn)用文中提出的模型在深海區(qū)域?qū)﹄妶?chǎng)偶極子進(jìn)行分析是可行的。

［1］威廉·H·海特，約翰·A·巴克.工程電磁場(chǎng)［M］.西安:西安交通大學(xué)出版社，2009.

［2］倪光正，楊仕友，錢秀英，等.工程電磁場(chǎng)數(shù)值計(jì)算［M］.北京:機(jī)械工業(yè)出版社，2006.

［3］Kasper C，April M G.Electroglalvanic Finite Element Analysis of Partially Protected Marine Structures［M］.CORROSION/82，Huaston，TX:NACE，1982:168.

［4］黃臨平，戴世坤.復(fù)雜條件下3D 電磁場(chǎng)有限元計(jì)算方法［J］.武漢:中國(guó)地質(zhì)大學(xué)學(xué)報(bào)，2002，17(4):373-377.

［5］劉勝道，龔沈光.并矢格林函數(shù)法求解海水中電偶極子電場(chǎng)［J］.電波科學(xué)學(xué)報(bào)，2002.

［6］LI Yong，LIN Pin rong，XU Bao li，et al.Computing the electromagnetic response of the time-harmonic horizontal electric dipole in a homog-enous half-space under ionosphere coupling［C］//China:Computing Techniques for Geophysical and Geochemical Exploration，2008.

［7］趙博，張洪亮.Ansoft12 在工程電磁場(chǎng)中的應(yīng)用［M］.北京:中國(guó)水利水電出版社，2011.