何　芳，王向軍，高俊吉

(海軍工程大學(xué)電氣工程學(xué)院,湖北 武漢　430033)

?

基于邊界元的艦船水下靜態(tài)電場(chǎng)推算方法

何芳，王向軍，高俊吉

(海軍工程大學(xué)電氣工程學(xué)院,湖北 武漢430033)

0引言

在艦船靜態(tài)電場(chǎng)建模方面，目前主要有邊界元法(BEM)[1-2]、有限元法(FEM)[3-4]、等效電偶極子法[5-6]，艦船電場(chǎng)的深度換算是電場(chǎng)建模的重要研究?jī)?nèi)容。艦船靜態(tài)電場(chǎng)的深度換算目前國(guó)外未見有報(bào)道，國(guó)內(nèi)可見報(bào)道很少，可查資料顯示只有海軍工程大學(xué)在該方面做了初步工作，主要采用點(diǎn)電源法[7-8]、電偶極子法[9-11]、平面迭代法[12]進(jìn)行深度換算。文獻(xiàn)[8]采用點(diǎn)電源產(chǎn)生的靜電場(chǎng)模擬船舶，利用泰勒展開式的收斂半徑對(duì)場(chǎng)源的個(gè)數(shù)和位置進(jìn)行了準(zhǔn)確的計(jì)算，建立了基于點(diǎn)電源的靜態(tài)模型，再反演得到模擬源的模型，最后利用此模型計(jì)算目標(biāo)深度上的靜電場(chǎng)。此方法對(duì)建模要求很高，模型的準(zhǔn)確度對(duì)換算結(jié)果影響較大，計(jì)算過程較為復(fù)雜。文獻(xiàn)[9-11]主要是基于電偶極子進(jìn)行無限大平面的深度換算，解決了電場(chǎng)某分量向下進(jìn)行較大深度自換算和互換算的問題，文獻(xiàn)[12]利用微分原理，根據(jù)靜態(tài)電場(chǎng)各分量與其垂直分量的偏導(dǎo)數(shù)關(guān)系，建立各無限大平面之間的數(shù)學(xué)關(guān)系，采用迭代算法解決海水中靜態(tài)電場(chǎng)向較淺深度換算問題。但是基于無限大平面理論的深度換算在轉(zhuǎn)化為數(shù)值積分后，對(duì)換算深度、測(cè)量點(diǎn)距及測(cè)量平面尺寸都有要求，若換算深度小于測(cè)量點(diǎn)距或測(cè)量平面尺寸接近換算深度時(shí)都會(huì)產(chǎn)生較大誤差，同時(shí)在工程實(shí)際應(yīng)用中，無法在無限大平面上測(cè)量，只能在有限平面內(nèi)取值，造成了換算精度的下降，不利用指導(dǎo)工程實(shí)踐。本文針對(duì)此問題，提出了基于邊界元的艦船水下靜態(tài)電場(chǎng)推算方法。

1水下靜態(tài)電場(chǎng)的邊界元模型

建立水下靜態(tài)電場(chǎng)邊界元模型的思路是：先建立半無限域任意形狀的閉合曲面邊界面積分模型；然后通過對(duì)邊界面進(jìn)行正方形面元剖分，將邊界積分方程離散為代數(shù)方程組，解出邊界源點(diǎn)近似解，并對(duì)邊界奇異9+++點(diǎn)進(jìn)行數(shù)值處理；最后根據(jù)域內(nèi)場(chǎng)點(diǎn)與邊界源點(diǎn)近似解的關(guān)系式，求解半無限域內(nèi)的靜態(tài)電場(chǎng)分布。

1.1任意形狀的閉合曲面邊界面積分模型

艦船靜態(tài)電場(chǎng)是由海水中的穩(wěn)恒電荷分布產(chǎn)生的，基于媒質(zhì)內(nèi)部均勻且各向同性的假設(shè),穩(wěn)恒電荷只分布在空氣和海水、海水和艦船以及海水和海床的分界面上, 海水內(nèi)部無凈電荷存在，電荷面密度與分界面兩邊媒質(zhì)的電容率有關(guān)。若在近場(chǎng)區(qū)取一任意形狀閉合曲面，海水無限遠(yuǎn)處邊界面視為S∞,S0-S∞構(gòu)成一無限域V,半無限域V的內(nèi)邊界面S0的外法線方向指向無限域內(nèi)部，如圖1所示。半無限域V內(nèi)的靜態(tài)電場(chǎng)則滿足Laplace方程。

(1)

設(shè)邊界面S0上有一源點(diǎn)M0,M為無限域V中的任意一場(chǎng)點(diǎn)，r為源點(diǎn)到場(chǎng)點(diǎn)的距離。半無限域內(nèi)場(chǎng)值與內(nèi)邊界面S0上函數(shù)值關(guān)系，可以通過Green第三公式[13]推導(dǎo)得出：

(2)

圖1　艦船水下電場(chǎng)半無限區(qū)域平面圖Fig.1　Schematic of vessel electric field in the semiinfinite domain

1.2邊界積分方程離散及奇異點(diǎn)處理

設(shè)邊界面S0上的有限測(cè)量點(diǎn)為M0，根據(jù)Green第二公式[13]推導(dǎo)出邊界積分公式為：

2πu(M0)=

(3)

(4)

式中，i=1,2,…,k。ri,j為面積元dSi的中心與面積元dSj的中心之間的距離。

(5)

當(dāng)i=j時(shí) ，場(chǎng)點(diǎn)和源點(diǎn)在同一積分面元上，邊界面積分出現(xiàn)奇異點(diǎn)。用正方形小面元(面元邊長(zhǎng)為a)，正方形面元坐標(biāo)設(shè)置如圖2所示，計(jì)算得

(6)

圖2　正方形面元的坐標(biāo)設(shè)置Fig.2　setting of coordinate in a element

用矩陣方程描述即為Ak×kU1×k=Bk×k,系數(shù)矩陣A，B均為k×k的滿秩方陣，U1×k即邊界面上的整個(gè)電場(chǎng)分布。為了便于說明，在笛卡兒坐標(biāo)系下對(duì)代數(shù)方程組系數(shù)A，B進(jìn)行確定。以平行于潛艇縱向方向由艦尾向艦首為X軸正方向，以平行于潛艇橫向方向由左舷向右舷為Y軸正方向，以垂直于水面指向地心為Z軸正方向。以潛艇為中心取一長(zhǎng)方體表面作為邊界閉合面，六個(gè)面上的法線方向都指向長(zhǎng)方體內(nèi)部，如圖3所示。

圖3　長(zhǎng)方體邊界面及坐標(biāo)系Fig.3　Boundary surface and its coordinate system

ri,j為面積元dSi的中心與面積元dSj的中心之間的距離：

(7)

只要邊界面的形狀和具體剖分方法確定了，矩陣方程Ak×k=Bk×k的系數(shù)矩陣A，B就能惟一確定。解此方程組,可以直接計(jì)算出U1×k，得到邊界面上各面元中心處的電勢(shì)u，即邊界面上源點(diǎn)的近似解。

1.3基于邊界元法的水下靜態(tài)電場(chǎng)分布

將邊界面上各面元中心處的電勢(shì)代入式(2)，得到半無限場(chǎng)域內(nèi)電場(chǎng)強(qiáng)度的函數(shù)解析解。

(8)

(9)

2基于邊界元的艦船水下靜態(tài)電場(chǎng)推算方法

為了接近實(shí)際情況(方便潛艇通過)，不考慮邊界面S1，S2兩個(gè)側(cè)面的電場(chǎng)分布，驗(yàn)證采用長(zhǎng)方體四個(gè)面近似閉合曲面的水下靜態(tài)電場(chǎng)邊界元模型是否準(zhǔn)確可靠，采用MATLAB進(jìn)行仿真。將一單位電偶極子放置于水深5 m處，設(shè)計(jì)邊界面為一長(zhǎng)方體外表面，如圖3所示，長(zhǎng)寬高分別為x∈(-6,6)m,y∈(-1.5,1.5)m,z∈(-1.5,1.5)m。忽略長(zhǎng)方體S1，S2兩個(gè)側(cè)面，將其余四個(gè)側(cè)面每個(gè)面都剖分為648個(gè)正方形面元，共計(jì)剖分2 592個(gè)面元，每個(gè)面元的中心點(diǎn)電場(chǎng)強(qiáng)度作為整個(gè)面元電場(chǎng)強(qiáng)度的平均值。 應(yīng)用式(9)，換算得到長(zhǎng)方體表面積外的半無限區(qū)域內(nèi)任意位置的電場(chǎng)三分量，以此作為水下靜態(tài)電場(chǎng)邊界元模型。

2.1利用邊界元模型外推直線處電場(chǎng)

利用水下靜態(tài)電場(chǎng)邊界元模型外推直線z=20 m,y=0,x∈(-50,50) m處的電場(chǎng)三分量，以及電偶極子在相同位置直接產(chǎn)生的電場(chǎng)三分量對(duì)照情況如圖4所示，與文獻(xiàn)[5]的電場(chǎng)分布解析表達(dá)式仿真結(jié)果(理論值)完全吻合。

2.2利用邊界元模型外推平面處電場(chǎng)

利用水下靜態(tài)電場(chǎng)邊界元模型外推平面z=20 m,y∈(-50,50) m,x∈(-50,50) m處的電場(chǎng)三分量，如圖5所示。仿真結(jié)果與文獻(xiàn)[9]采用無限大平面法換算得到的相同平面處電場(chǎng)三分量在數(shù)量級(jí)上是完全吻合的。

仿真結(jié)果證明了采用長(zhǎng)方體四個(gè)面近似閉合曲面的邊界元模型是可行的，且驗(yàn)證了基于該邊界元模型的水下靜態(tài)電場(chǎng)推算方法的正確性和可靠性。

圖4　電場(chǎng)三分量仿真推算值和理論值對(duì)比圖Fig.4　The results of the calculated and the theoretical for the components of electric field intensity

圖5　邊界元換算平面電場(chǎng)分布圖Fig.5　Result of BEM on plane

3實(shí)驗(yàn)驗(yàn)證

為了驗(yàn)證基于邊界元水下靜態(tài)電場(chǎng)推算方法的有效性及換算精度，在實(shí)驗(yàn)室采用潛艇模型進(jìn)行實(shí)驗(yàn)驗(yàn)證。

3.1實(shí)驗(yàn)條件和實(shí)驗(yàn)步驟

實(shí)驗(yàn)硬件設(shè)施包括無磁性實(shí)驗(yàn)水池、實(shí)驗(yàn)潛艇模型、電場(chǎng)三分量傳感器及矩形測(cè)量架。無磁性實(shí)驗(yàn)水池，尺寸為8m×5m×1.5m,水深0.8m，將工業(yè)用鹽倒入池水中，充分溶解并混合均勻，模擬淡鹽水，測(cè)得所配制的海水電導(dǎo)率為3.96S/m。潛艇模型長(zhǎng)1.5m,寬0.13m，高0.15m，按實(shí)艇縮小比例1∶100制得。艇模船殼鐵制，螺旋槳青銅制，船模的不同金屬結(jié)構(gòu)在模擬淡鹽水中因電化學(xué)作用產(chǎn)生腐蝕電流，為防止船殼腐蝕，艇模上裝載了陰極保護(hù)系統(tǒng)(ICCP)，而防腐系統(tǒng)的開啟也會(huì)在模擬淡鹽水中產(chǎn)生電流。這些腐蝕和防腐電流在模擬海水(池水)中形成了靜電場(chǎng)。

矩形測(cè)量架寬0.4m,高0.4m。完全置于水中固定，測(cè)量架每條邊上均勻布放3個(gè)電場(chǎng)三分量傳感器。

實(shí)驗(yàn)步驟如下：

1)將艇模移至距測(cè)量架較遠(yuǎn)處，電場(chǎng)三分量傳感器采集背景電場(chǎng)。

2)以測(cè)量架為中心建立笛卡兒坐標(biāo)系，潛艇沿X軸勻速穿過測(cè)量架，如圖6所示(潛艇由紙內(nèi)向紙面運(yùn)動(dòng))，電場(chǎng)三分量傳感器可以采集到艇模在各個(gè)位置時(shí)的電場(chǎng)值，即可得到一個(gè)除前后兩個(gè)端面外的矩形包絡(luò)面。用矩形包絡(luò)面上各測(cè)量點(diǎn)處的電場(chǎng)減去背景電場(chǎng)即可得到艇模在矩形包絡(luò)面上的電場(chǎng)分布。

圖6　潛艇與傳感器位置示意圖Fig.6　Schematic of position of submarine and sensor

3)根據(jù)所建立的靜態(tài)電場(chǎng)邊界元模型推算距潛艇左舷3m,水深0.2m,即直線 z=0.2m,y=3m,x∈(-8,8)m處的電場(chǎng)。

圖7　電場(chǎng)強(qiáng)度測(cè)量值與邊界元推算值的比較Fig.7　The comparision of electric field intensity betweenthe measured and the computed

3.2實(shí)測(cè)數(shù)據(jù)結(jié)果分析

根據(jù)等效偶極子理論[5-6]，實(shí)驗(yàn)中船模可以看做電場(chǎng)源點(diǎn)，z=0.2m,y=3m,x∈(-8,8)m直線處各點(diǎn)作為場(chǎng)點(diǎn)。以船模中心為對(duì)稱軸，在偏離源點(diǎn)較近處電場(chǎng)模值出現(xiàn)大幅度衰減，偏離源點(diǎn)較遠(yuǎn)處衰減變?nèi)?。在該?shí)驗(yàn)建立的笛卡兒坐標(biāo)系中，船模中心坐標(biāo)(0，0，0)，場(chǎng)點(diǎn)坐標(biāo)(0，3，0.2)及(6，3，0.2)代入文獻(xiàn)[5]中靜電場(chǎng)分布解析公式驗(yàn)證，可以發(fā)現(xiàn)在船長(zhǎng)約4倍處，靜態(tài)電場(chǎng)強(qiáng)度模量衰減約75%，4倍以內(nèi)衰減較快，4倍以外衰減變緩。該實(shí)驗(yàn)實(shí)測(cè)數(shù)據(jù)與文獻(xiàn)[5]中直流電偶極子等效理論公式在相同源點(diǎn)與場(chǎng)點(diǎn)處計(jì)算結(jié)果一致，證明該實(shí)驗(yàn)是有效的。

3.3邊界元模型換算數(shù)據(jù)分析

根據(jù)式(9)可知，電場(chǎng)強(qiáng)度與場(chǎng)點(diǎn)和源點(diǎn)之間的距離成反比。當(dāng)場(chǎng)點(diǎn)與源點(diǎn)距離較近時(shí)，電場(chǎng)強(qiáng)度衰減越快；場(chǎng)點(diǎn)與源點(diǎn)距離較遠(yuǎn)時(shí)，電場(chǎng)強(qiáng)度衰減越慢。定性上來看，靜態(tài)電場(chǎng)邊界元模型推算的電場(chǎng)強(qiáng)度模值是符合靜電場(chǎng)在海水中的衰減規(guī)律的。根據(jù)邊界元模型換算數(shù)據(jù)顯示，在船長(zhǎng)約4倍處，靜態(tài)電場(chǎng)強(qiáng)度模量衰減約75%，與實(shí)驗(yàn)實(shí)測(cè)數(shù)據(jù)衰減率接近，可以把船長(zhǎng)4倍處看做電場(chǎng)強(qiáng)度模值衰減特征點(diǎn)。

4結(jié)論

本文提出了基于邊界元的艦船水下靜態(tài)電場(chǎng)推算方法。該方法不局限于無限大平面，可對(duì)半無限域任意形狀的閉合曲面建立邊界面積分模型，并進(jìn)一步采用長(zhǎng)方體四個(gè)面近似閉合曲面建立水下靜態(tài)電場(chǎng)邊界元模型，基于該邊界元模型進(jìn)行水下靜態(tài)電場(chǎng)的推算。實(shí)驗(yàn)驗(yàn)證和結(jié)果分析證明了基于邊界元的水下靜態(tài)電場(chǎng)推算方法的正確性和可靠性，其外推電場(chǎng)值與實(shí)測(cè)電場(chǎng)值相對(duì)誤差為5%左右，能達(dá)到較高的換算精度，相比于無限大平面法的深度換算，更適宜指導(dǎo)工程實(shí)踐。

基于邊界元的艦船水下靜態(tài)電場(chǎng)推算方法可以在以下方面做進(jìn)一步改進(jìn)：1)基于實(shí)際測(cè)量的原因，本文忽略了包絡(luò)面的前、后兩個(gè)端面，來近似封閉面，由于前后兩端面面積較小，仿真結(jié)果表明誤差不大，但仍會(huì)產(chǎn)生一定的計(jì)算誤差。2)在對(duì)邊界面進(jìn)行剖分時(shí)，將每個(gè)面元中心點(diǎn)上的值作為該面元上的平均值，這會(huì)產(chǎn)生離散誤差，因此小面元剖分得越細(xì)，面元越多，面元的面積越小，邊界元推算的精度越高。3)當(dāng)場(chǎng)點(diǎn)趨于源點(diǎn)時(shí)，邊界面積分涉及到對(duì)奇異點(diǎn)的數(shù)值處理，邊界面積分精度很大程度上依賴于奇異積分的計(jì)算精度，可以在對(duì)奇異積分處理方面對(duì)文中的邊界元模型進(jìn)行改進(jìn)。

參考文獻(xiàn):

[1]ZAMANING.Boundaryelementsimulationofthecathodicprotectionsysteminaprototypeship[J].AppliedMathematicsandComputation, 1988, 26:119-123.

[2]IWATAM,HUANGY,FUJIMOTOY.ApplicationofBEMtodesignoftheimpressedcurrentcathodicprotectionsystemforshiphull[J].JoftheSocietyofNavalArchitectsofJapan, 1992, 171: 377-380.

[3]DOIGP,FLFWITTPEJ.Afinitedifferencenumericalanalysisofgalvaniccorrosion[J].Electrochem,1979,126(12):2057-2063.

[4]KASPERRG,APRILMG.Electrogalvanicfiniteelementanalysisofpartiallyprotectedmarinestructures[J].Corrosion, 1983, 39(5): 181-183.

[5]劉勝道.艦船水下電場(chǎng)的測(cè)試技術(shù)與電偶極子模型研究[D].武漢：海軍工程大學(xué), 2002：43-45.

[6] 劉文寶，王向軍，嵇斗.基于電偶極子模型的艦船靜電場(chǎng)深度換算[J].空軍雷達(dá)學(xué)院學(xué)報(bào)，2010(6)：435-438.

[7] 劉忠樂，龔沈光.海水中穩(wěn)恒電流電場(chǎng)的點(diǎn)電極計(jì)算模型[J].海軍工程大學(xué)學(xué)報(bào)，2004，16(1)：35-39.

[8] 姜潤(rùn)翔，胡英娣，龔沈光.基于點(diǎn)電源的船舶靜態(tài)電場(chǎng)深度換算方法研究[J].電波科學(xué)學(xué)報(bào)，2014，29(4)：685-693.

[9] 陳聰，李定國(guó)，龔沈光.基于拉氏方程的艦船靜態(tài)電場(chǎng)深度換算[J].電子學(xué)報(bào)，2010，38(9)：2025-2030.

[10]李定國(guó)，陳聰，劉華波.艦船水下靜態(tài)電場(chǎng)深度換算的格林函數(shù)法[J].海軍工程大學(xué)學(xué)報(bào)，2012，24(3)：1-5.

[11] 陳聰，龔沈光，李定國(guó).艦船靜態(tài)電場(chǎng)深度換算方法[J].哈爾濱工程大學(xué)學(xué)報(bào)，2009，30(6)：719-722.

[12]胡英娣，龔沈光，閆永貴.一種新的船舶靜態(tài)電場(chǎng)深度換算方法[J].海軍工程大學(xué)學(xué)報(bào)，2013，10(25)：17-20.

[13]姚端正，梁家寶.數(shù)學(xué)物理方法[M].武漢：武漢大學(xué)出版社，2010:250-252.

摘要:針對(duì)基于無限大平面理論的深度換算在轉(zhuǎn)化為數(shù)值積分后，對(duì)換算深度、測(cè)量點(diǎn)距及測(cè)量平面尺寸都有要求，若換算深度小于測(cè)量點(diǎn)距或測(cè)量平面尺寸接近換算深度時(shí)都會(huì)產(chǎn)生較大誤差等問題，提出了基于邊界元的艦船水下靜態(tài)電場(chǎng)推算方法。該方法不局限于無限大平面，可對(duì)半無限域任意形狀的閉合曲面建立邊界面積分模型，進(jìn)一步采用長(zhǎng)方體四個(gè)面近似閉合曲面建立水下靜態(tài)電場(chǎng)邊界元模型，并利用該邊界元模型進(jìn)行了水下靜態(tài)電場(chǎng)的推算。實(shí)驗(yàn)驗(yàn)證和結(jié)果分析證明了此方法的有效性，能達(dá)到較高的換算精度。

關(guān)鍵詞:邊界元；靜態(tài)電場(chǎng)；半無限區(qū)域

Extrapolation of Vessel Static Electric Field Based on Boundary ElementHE Fang，WANG Xiangjun，GAO Junji

(Electrical Engineering College，Naval University of Engineering，Wuhan 430033，China)

Abstract:When the depth is less than?measuring distance, or the size of measuring plane is close to the depth, extrapolation based on the infinite plane will have a greater error. If extrapolation is based on the finite plane, a decrease in accuracy is caused. To solve this problem, a method of depth extrapolation of static electric field was proposed based on boundary element method in the semi infinite domain. The method was not limited to the infinite plane, but also could establish the boundary integral model for arbitrary closed surface in semi infinite domain. Static electric field boundary element model under water was built up based on four surface of rectangular approximated the closed surface, and the underwater static electric field could be calculated using this boundary element model. The submarine model experiment testified that this method was of higher precision and stronger operability.

Key words:boundary element; static electric field; semi infinite domain

中圖分類號(hào):TM46

文獻(xiàn)標(biāo)志碼:A

文章編號(hào)：1008-1194(2015)06-0055-06

作者簡(jiǎn)介:何芳(1981—), 女,湖北荊門人，講師, 博士研究生，研究方向：電場(chǎng)防護(hù)與隱身技術(shù)。E-mail:wskay1103@163.com。

基金項(xiàng)目:國(guó)家自然科學(xué)基金項(xiàng)目資助(41476153)

*收稿日期:2015-03-27