黨偉+孫雨

摘要：在運(yùn)用邊界元方法計(jì)算開(kāi)域場(chǎng)中的電場(chǎng)時(shí)，為提高傳統(tǒng)一階邊界元法的精度問(wèn)題，提出NURBS曲面四邊形邊界元法。此方法是對(duì)計(jì)算模型進(jìn)行二階四邊形剖分，采用剖分信息構(gòu)造出NURBS曲面參數(shù)方程，然后用面積比值法定義曲面上頂點(diǎn)的形狀函數(shù)，曲面四邊形的四個(gè)頂點(diǎn)作為求解點(diǎn)，計(jì)算NURBS曲面四邊形單元上的積分值。

關(guān)鍵詞：NURBS；電場(chǎng)強(qiáng)度；提高精度；電場(chǎng)計(jì)算；邊界元法 文獻(xiàn)標(biāo)識(shí)碼：A

中圖分類號(hào)：TM151 文章編號(hào)：1009-2374（2016）29-0018-04 DOI：10.13535/j.cnki.11-4406/n.2016.29.008

邊界元法在求解電磁數(shù)值計(jì)算中，是利用離散技術(shù)將求解的邊界面離散為邊界單元，將邊界問(wèn)題等價(jià)轉(zhuǎn)化為邊界積分方程問(wèn)題，同時(shí)把邊界積分方程離散為代數(shù)方程組，最后利用數(shù)值方法計(jì)算方程組，得到邊界積分方程的解。

本文提出NURBS曲面四邊形邊界元方法。利用二階剖分單元的節(jié)點(diǎn)信息，構(gòu)造出NURBS曲面四邊形單元，以曲面四邊形單元的頂點(diǎn)為求解點(diǎn)，計(jì)算NURBS曲面四邊形單元上的積分，最后形成NURBS曲面四邊形邊界元方法求解三維電場(chǎng)邊值問(wèn)題。

1 NURBS曲面邊界元的基本原理

1.1 有理B樣條基函數(shù)

NURBS（Non-Uniform Rational B-Spline）方法，是在B樣條方法基礎(chǔ)上發(fā)展起來(lái)的，為了能進(jìn)一步精確地描述有理數(shù)，推廣得到有理多項(xiàng)式基函數(shù)。將B樣條基函數(shù)推廣到有理的基函數(shù)，就得到了NURBS方法的基函數(shù)即有理B樣條基函數(shù)，而且NURBS曲面在中心投影下得到B樣條曲面，所以研究NURBS方法就必須了解B樣條方法基礎(chǔ)。

B樣條函數(shù)的概念最初是由Schoenberg于1946年提出的。由de Boor-Cox公式得出提出B樣條的遞推公式為：

式中：稱為相應(yīng)于序列的次，即階B樣條基函數(shù)。其中，下標(biāo)表示序號(hào)，下標(biāo)表示次數(shù)?？梢钥闯鱿胍玫降趥€(gè)次B樣條基函數(shù)，就得共個(gè)節(jié)點(diǎn)序列，而且的支撐區(qū)間為。由該定義出發(fā)，可以導(dǎo)出任意階次的B樣條基函數(shù)。文章中是采用高斯積分方法求解數(shù)值積分，由于高斯積分法對(duì)應(yīng)的積分區(qū)間為。而在區(qū)間上B樣條基函數(shù)為均勻的基函數(shù)，根據(jù)B樣條基函數(shù)的性質(zhì)，次B樣條基函數(shù)只有個(gè)基函數(shù)在上非零，而且在上都為多項(xiàng)式函數(shù)。取，則由式（1）可得均勻2次B樣條基函數(shù)表達(dá)式為：

有理B樣條基函數(shù)作為對(duì)B樣條基函數(shù)的推廣，一個(gè)基函數(shù)變?yōu)橛欣矸质绞紫刃枰鄳?yīng)的一個(gè)多項(xiàng)式作為分母，然而考慮到相同的分母將會(huì)給有理分式之間的計(jì)算帶來(lái)簡(jiǎn)便。再注意到B樣條基函數(shù)具有的“權(quán)”性，即在區(qū)間的非負(fù)性和單位分解性，因此可以采取下面的方法：給每一個(gè)B樣條基函數(shù)乘以一個(gè)“權(quán)因子”，再用它們的“和”作為統(tǒng)一分母進(jìn)行“平均”，這樣就得到了次有理B樣條基函數(shù)為：

1.2 NURBS曲面

NURBS曲面是非有理張量積型性B樣條曲面的有理推廣，類似于NURBS曲線定義，一張次NURBS曲面有理式可以表示為：

其中是雙變量的有理基函數(shù)：

式中：為控制點(diǎn)；為與控制點(diǎn)相應(yīng)的權(quán)因子；和分別為參數(shù)向次和參數(shù)向次的定義在和規(guī)范有理B樣條基函數(shù)。

根據(jù)式（3），當(dāng)時(shí)，得到廣泛應(yīng)用的雙二次NURBS曲面的數(shù)學(xué)模型為：

上式中是由式（2）推導(dǎo)出，為的轉(zhuǎn)置矩陣。

根據(jù)式（7），要得出NURBS曲面的參數(shù)方程，必須就算出控制點(diǎn)P11，P12，…，P33。根據(jù)有理式表示的NURBS曲面參數(shù)方程還可以表示為，其中為曲面數(shù)據(jù)點(diǎn)列，維數(shù)為；為有理基函數(shù)，維數(shù)為，為方陣，則逆矩陣存在；是控制點(diǎn)，維數(shù)為。在已知曲面上的數(shù)據(jù)點(diǎn)列，假定權(quán)因子時(shí)，經(jīng)過(guò)矩陣求逆變換可得到雙二次NURBS曲面控制點(diǎn)矩陣求逆變換可得。

1.3 四邊形單元上的形狀函數(shù)

本文是以四邊形單元的四個(gè)頂點(diǎn)為求解點(diǎn)，運(yùn)用面積比值法構(gòu)造曲面四邊形頂點(diǎn)上的形狀函數(shù)。曲面單元如圖1所示，曲面單元上任意一點(diǎn)處的函數(shù)值由4個(gè)頂點(diǎn)節(jié)點(diǎn)i、j、k、l上的函數(shù)值和點(diǎn)相對(duì)于4個(gè)頂點(diǎn)節(jié)點(diǎn)的形狀函數(shù)插值得到。在圖1中，點(diǎn)處的兩條參數(shù)等值線將曲面單元分成、、、共4部分，其中每一部分的面積都可以通過(guò)曲面參數(shù)方程得到。設(shè)曲面單元總的面積為。則點(diǎn)處關(guān)于4個(gè)頂點(diǎn)的形狀函數(shù)定義為：

1.4 邊界元法

三維靜電場(chǎng)光滑邊界的積分方程為：

式中：為邊界上的場(chǎng)點(diǎn)；為邊界上的源點(diǎn)；為場(chǎng)點(diǎn)電位；為源點(diǎn)電位；為源點(diǎn)到場(chǎng)點(diǎn)的距離，；為邊界外法線方向；S為求解區(qū)域的邊界。

2 NURBS曲面邊界元算例分析

2.1 算例一

如圖2是高為2、半徑為1的四分之一圓柱面。運(yùn)用面積計(jì)算可知其表面積為。將其剖分一個(gè)單元，得出節(jié)點(diǎn)信息后，導(dǎo)入一階平面四邊形邊界元程序計(jì)算其表面積為2.8284，相對(duì)誤差為9.97%。將節(jié)點(diǎn)信息導(dǎo)入NURBS曲面邊界元程序計(jì)算出四分之一圓柱面的表面積為3.1645，相對(duì)誤差為0.073%。可以得出NURBS曲面邊界元程序相比一階平面四邊形邊界元程序精度提高了9.24%，所以NURBS曲面更接近實(shí)際的四分之一圓柱面，NURBS曲面邊界元比一階平面四邊形邊界元計(jì)算更加精確。

2.2 算例二

半徑為1的導(dǎo)體球面，給導(dǎo)體球面加上100V電壓，求解球面上的電場(chǎng)強(qiáng)度。用ANSYS建立球面模型，將其自由剖分54個(gè)單元，導(dǎo)出單元及節(jié)點(diǎn)信息代入NURBS曲面邊界元程序計(jì)算各節(jié)點(diǎn)上的電場(chǎng)強(qiáng)度，結(jié)果云圖如圖3（a）。建立模型后進(jìn)行一階單元剖分，自由剖分為54個(gè)單元，代入一階平面四邊形邊界元程序計(jì)算得出結(jié)果為圖3（b）。將兩種計(jì)算沿球面電場(chǎng)強(qiáng)度的結(jié)果在同一坐標(biāo)系中呈現(xiàn)為圖3（c）。

由電磁場(chǎng)知識(shí)得知半徑為1m的導(dǎo)體球面加上100V電壓，則球面上的電場(chǎng)強(qiáng)度為100V/m。NURBS曲面邊界元計(jì)算結(jié)果最大相對(duì)誤差為1.46%，一階平面四邊形邊界元計(jì)算結(jié)果最大誤差為12.94%。從兩種計(jì)算得出沿球面上電場(chǎng)強(qiáng)度分布，可以得出NURBS曲面邊界元計(jì)算結(jié)果基本上落在電場(chǎng)強(qiáng)度100V/m左右，節(jié)點(diǎn)分布相對(duì)規(guī)律。一階平面四邊形邊界元計(jì)算結(jié)果距離電場(chǎng)強(qiáng)度100V/m偏差加大，節(jié)點(diǎn)分布相對(duì)零亂。所以可以得出在相同的剖分單元下，NURBS曲面單元擬合出的曲面更接近實(shí)際的球面，計(jì)算結(jié)果相對(duì)一階平面也更加精確。

2.3 算例三

陶瓷絕緣子型號(hào)XP-120，絕緣子結(jié)構(gòu)高度為146mm，建立該絕緣子模型后對(duì)邊界面采用二階四邊形單元規(guī)則剖分240個(gè)單元，如圖4（a）。高壓端金具施加kV電壓，低壓端金具電壓為0V。傘群和芯棒認(rèn)為是一種介質(zhì)，相對(duì)介電常數(shù)為6.0，分別用NURBS曲面邊界元法和一階平面四邊形邊界元法計(jì)算絕緣子的電場(chǎng)強(qiáng)度。為了驗(yàn)證邊界元算法計(jì)算結(jié)果的正確性，根據(jù)軸對(duì)稱可以用二維軸對(duì)稱高精度有限元法（無(wú)限遠(yuǎn)處為零電位）計(jì)算絕緣子電場(chǎng)強(qiáng)度進(jìn)行結(jié)果對(duì)比分析驗(yàn)證，單位為kV/mm。

由圖4從數(shù)值可以看出，NURBS曲面邊界元法計(jì)算的絕緣子上最大電場(chǎng)強(qiáng)度為13.95kV/mm，高精度有限元計(jì)算結(jié)果為13.134kV/mm，一階平面四邊形邊界元法計(jì)算結(jié)果為17.675kV/mm。對(duì)比圖4（b）、圖4（c）、圖4（d）可知，NURBS曲面邊界元法與高精度有限元算法計(jì)算結(jié)果基本一致；一階平面四邊形邊界元法計(jì)算結(jié)果相對(duì)誤差為34.5%左右，可以得出NURBS曲面邊界元法比一階平面四邊形邊界元法明顯的精度提高。邊界元算法只需要對(duì)求解模型的邊界進(jìn)行剖分，求解的規(guī)模明顯小一些，這也是邊界元法相比有限元法的優(yōu)勢(shì)。

3 結(jié)語(yǔ)

相比于傳統(tǒng)的邊界元法，提出了NURBS曲面四邊形邊界元方法，該方法在不增加節(jié)點(diǎn)的情況下提高邊界元方法計(jì)算的精度。通過(guò)算例表明，該方法比傳統(tǒng)一階平面單元更逼近實(shí)際邊界，既提高了精度，又沒(méi)有增加節(jié)點(diǎn)數(shù)，達(dá)到了比較理想的效果。

參考文獻(xiàn)

[1] 姚振漢，王海濤.邊界元法[M].北京：高等教育出版社，2010.

[2] 王仁宏，李崇軍，朱春鋼.計(jì)算幾何教程[M].北京：科學(xué)教育出版社，2008.

[3] 施法中.計(jì)算機(jī)輔助幾何設(shè)計(jì)與非均勻有理B樣條[M].北京：高等教育出版社，2001.

[4] 李亞莎，王澤忠.基于圓環(huán)坐標(biāo)系的三維靜電場(chǎng)曲面三角形邊界元法[J].電工技術(shù)學(xué)報(bào)，2006，21（9）.

作者簡(jiǎn)介：黨偉（1975-），男，湖北十堰人，國(guó)網(wǎng)湖北省電力公司十堰供電公司高級(jí)工，助理工程師。

（責(zé)任編輯：黃銀芳）