王韶霈，宋開(kāi)宏，謝國(guó)大，黃 帥

（安徽大學(xué) 電子信息工程學(xué)院，安徽 合肥 230601）

0 引 言

在電磁場(chǎng)工程問(wèn)題中，靜電場(chǎng)和時(shí)諧場(chǎng)的場(chǎng)分布是電磁場(chǎng)工程設(shè)計(jì)首先要解決的問(wèn)題，因?yàn)樗苯雨P(guān)系到整個(gè)系統(tǒng)的性能是否滿足實(shí)際要求。隨著對(duì)靜電現(xiàn)象研究的日益深入，對(duì)靜電場(chǎng)中場(chǎng)強(qiáng)和電位分布的了解也更加迫切。隨著通信業(yè)、航空航天業(yè)的發(fā)展，對(duì)高空中電離層的研究有著重要的意義。電離層作為一種傳播介質(zhì)會(huì)使電波損失部分能量，對(duì)通信造成一定影響。就以上問(wèn)題，需要在地面建立帶有電極的真空罐［1］對(duì)其進(jìn)行模擬仿真，觀察研究其特性。對(duì)真空罐的研究也要先從靜電場(chǎng)做起，對(duì)靜電場(chǎng)進(jìn)行研究就需要通過(guò)描繪出其電位分布來(lái)具體分析。

在電磁場(chǎng)數(shù)值分析計(jì)算方法中，有限差分法一直是一種重要的數(shù)值方法。靜電場(chǎng)中電位函數(shù)滿足拉普拉斯方程，通過(guò)給定的邊界條件求解拉普拉斯方程，最終得到靜電場(chǎng)內(nèi)的電位分布圖以及電場(chǎng)圖。在對(duì)待二維靜電場(chǎng)問(wèn)題中運(yùn)用有限差分法求拉普拉斯方程數(shù)值解時(shí)最為經(jīng)典是運(yùn)用五點(diǎn)差分格式。對(duì)待三維靜電場(chǎng)問(wèn)題時(shí)會(huì)應(yīng)用到七點(diǎn)差分格式，若所求區(qū)域?yàn)檩S對(duì)稱，七點(diǎn)差分格式可以轉(zhuǎn)化為五點(diǎn)差分格式。這種利用離散節(jié)點(diǎn)與周圍四個(gè)節(jié)點(diǎn)的關(guān)系迭代求解電位值的方法很適合在EXCEL上完成。本文從五點(diǎn)差分格式出發(fā)，運(yùn)用EXCEL計(jì)算真空罐里的靜電場(chǎng)的電位數(shù)值，聯(lián)合MATLAB的圖形功能，繪制靜電場(chǎng)等位線，形象直觀掌握真空罐中的場(chǎng)分布，為工程設(shè)計(jì)奠定基礎(chǔ)。

1 理論分析

柱坐標(biāo)下的拉普拉斯方程有以下形式：

設(shè)空間中任意一點(diǎn)為U（z0，r0，θ0），根據(jù)有限差分法［2］，取圓柱等間距網(wǎng)格高Δz＝h，半徑Δr＝k，角度Δθ＝τ。設(shè)：

圓柱坐標(biāo)系中三維拉普拉斯方程的七點(diǎn)差分公式［3］：

其中：

C0、C1、C3、C4、C5、C6是取決于步長(zhǎng)的已知系數(shù)，U0、U1、U2、U3、U4、U5、U6則是待求的未知數(shù)。

待求區(qū)域如果具有對(duì)稱性，這個(gè)區(qū)域內(nèi)各點(diǎn)的電位函數(shù)值也與角度θ無(wú)關(guān)，且介質(zhì)均勻，則式（2）有關(guān)對(duì)應(yīng)角度θ的項(xiàng)就可不用考慮，那么七點(diǎn)差分方程就可簡(jiǎn)化為五點(diǎn)差分方程：

而r0＝0處即對(duì)稱軸上，式（3）不適用，對(duì)稱軸上的點(diǎn)系數(shù)為：

從式（3）可看出，轉(zhuǎn)化后的公式計(jì)算僅在rz平面，就位于rz平面內(nèi)的軸對(duì)稱平面而言，它的場(chǎng)域離散化下的網(wǎng)格劃分與離散點(diǎn)分布方式與二維場(chǎng)中的討論相同。在軸對(duì)稱場(chǎng)中，網(wǎng)格在空間的延續(xù)分布，即是關(guān)于軸對(duì)稱的一系列回轉(zhuǎn)體。

2 實(shí)際工程問(wèn)題的計(jì)算機(jī)解法

2.1 具體問(wèn)題

在一個(gè)直徑為70 cm、長(zhǎng)度為180 cm、兩端為半球的真空罐內(nèi)，在其軸心中間對(duì)稱、間距為60 cm的位置放置兩個(gè)環(huán)形電極，電極上環(huán)形電極的直徑為4 cm，不考慮電極環(huán)的厚度尺寸，真空罐按接地處理，環(huán)狀電極施加電壓信號(hào)為10 V，計(jì)算真空罐內(nèi)部空間的電場(chǎng)分布。

2.2 問(wèn)題的解決

真空罐內(nèi)空間是具有對(duì)稱性的，且介質(zhì)均勻，根據(jù)以上理論分析可以將這個(gè)三維問(wèn)題轉(zhuǎn)化為二維問(wèn)題，只計(jì)算分析真空罐一個(gè)截面里的電位分布，根據(jù)（3）式，用高斯—賽德?tīng)柕ǎ?］計(jì)算，公式為：

2.2.1 選取計(jì)算工具

根據(jù)式（4）可以看出，在數(shù)值求解時(shí)通常需要編程計(jì)算機(jī)完成大規(guī)模遞推關(guān)系，必須求解大型多元方程組，計(jì)算量龐大，給應(yīng)用中造成很大的不方便。EXCEL表在數(shù)值計(jì)算中有很大的應(yīng)用潛能，公式的復(fù)制十分突出地顯示出它的遞推運(yùn)算優(yōu)勢(shì)，特別是“重復(fù)計(jì)算”功能可十分簡(jiǎn)便地處理一些復(fù)雜遞推關(guān)系的數(shù)值計(jì)算，一般應(yīng)用軟件無(wú)法比擬。根據(jù)五點(diǎn)差分格式以及計(jì)算過(guò)程中需要多步迭代的特點(diǎn)，EXCEL在計(jì)算上有很大的優(yōu)勢(shì)，選擇［5，6］其為計(jì)算工具。

2.2.2 網(wǎng)格剖分

根據(jù)已知條件，本文在EXCEL中取71行與181列，同時(shí)要對(duì)真空罐截面的邊緣進(jìn)行數(shù)字化。圓弧邊緣作以下處理：取EXCEL中兩行36列單元格，從A1單元格依次輸入0～35，再在B1單元格中輸公式“35＋（35＾2－（A1－35）＾2）＾（1／2）”，將B1單元格的公式復(fù)制到后面35個(gè)單元格中，再對(duì)通過(guò)公式得到的數(shù)據(jù)四舍五入取整，最終獲得以（35，35）為圓心，35為半徑圓的四分之一圓邊在EXCL的坐標(biāo)，然后通過(guò)對(duì)稱關(guān)系補(bǔ)齊真空罐截面另外的圓弧邊界。

2.2.3 計(jì)算結(jié)果

剖分完成后的一個(gè)單元格即代表求解區(qū)域的一個(gè)節(jié)點(diǎn)，而單元格里的值就是該節(jié)點(diǎn)的電勢(shì)，這時(shí)就需要確定邊界條件，將真空罐罐體上的單元格都輸入數(shù)值0。由于電極厚度不計(jì)，則其在截面上可用2個(gè)點(diǎn)來(lái)表示，把代表兩個(gè)電極位置的單元格找到，輸入數(shù)值10。真空罐以外的單元格統(tǒng)一輸入數(shù)值0。

確定了邊界條件，需要輸入公式進(jìn)行計(jì)算。根據(jù)式（4），則每列的C3、C4是不同的，取決于r0，需要逐行輸入公式從對(duì)稱軸I，r0＝0開(kāi)始，每向上或者向下一行r0的值加1。逐行輸入公式時(shí)，選中需要輸入公式的一行，再選中本行其中一個(gè)單元格輸入相應(yīng)的公式，然后再選中這個(gè)單元格，復(fù)制，按住Ctrl鍵拖拽選取包含本行里的其他單元格，不要選中已經(jīng)輸入邊界條件的單元格，最后粘貼。進(jìn)行迭代計(jì)算，打開(kāi)“工具”菜單下的“選項(xiàng)”，選擇“重復(fù)計(jì)算”，將“迭代計(jì)算”選中，設(shè)置“最多迭代次數(shù)”為3 000，最大誤差為0.001，點(diǎn)擊“確定”，通過(guò)一段時(shí)間的迭代計(jì)算，得出結(jié)果，這時(shí)可以很直觀地看出計(jì)算結(jié)果。

由于EXCEL繪圖功能有限，而MATLAB擁有強(qiáng)大的繪圖功能，能彌補(bǔ)繪圖上的這一缺點(diǎn)，所以在數(shù)據(jù)到圖形的處理工作上，本文將選用MATLAB軟件來(lái)進(jìn)行。打開(kāi)MATLAB，利用文件菜單下的“Import data...”將EXCEL中計(jì)算得到的最終數(shù)據(jù)導(dǎo)入MATLAB，將導(dǎo)入的數(shù)據(jù)名定義為A。再用contour函數(shù)繪制等位線，在命令窗口輸入：contour（A，11），繪制了11條等位線，每條等位線間隔為1 V，這樣就得到了真空罐一個(gè)離散截面的電場(chǎng)等位線圖。

2.3 結(jié)論驗(yàn)證

MATLAB是一種功能十分強(qiáng)大、運(yùn)算效率很高的專業(yè)計(jì)算機(jī)軟件，它已經(jīng)成為一種極其靈活的計(jì)算機(jī)體系，具有簡(jiǎn)單易學(xué)、短小高效的源代碼，豐富的內(nèi)部函數(shù)等優(yōu)點(diǎn)。本文選用MATLAB對(duì)真空罐進(jìn)行仿真，并對(duì)結(jié)論進(jìn)行驗(yàn)證。同樣以式（4）為計(jì)算的基本公式，編寫程序并畫(huà)圖，得到了真空罐里的電位仿真圖，與運(yùn)用EXCEL計(jì)算得到的最終電位圖作比較分析，結(jié)果一致。

3 結(jié) 論

本文將有限差分法運(yùn)用到解決實(shí)際工程問(wèn)題中，這種以有限差分法為核心理論的方法相當(dāng)于一種仿真模擬軟件。它在計(jì)算二維靜電場(chǎng)問(wèn)題上受到的條件限制不多，只要網(wǎng)格剖分足夠精細(xì)，都能計(jì)算不規(guī)則邊界的靜電場(chǎng)問(wèn)題。而在三維問(wèn)題上目前只能計(jì)算軸對(duì)稱場(chǎng)的工程問(wèn)題，所以在計(jì)算三維問(wèn)題時(shí)還是有很大的局限性。雖然如此，但對(duì)EXCEL的應(yīng)用仍然有著方便快捷、計(jì)算結(jié)果以及通過(guò)它得到的圖像直觀明了、不需要編寫程序的優(yōu)勢(shì)，在解決更多工程問(wèn)題上還有很大的潛能等待開(kāi)發(fā)。

用有限差分法代替偏微分求解電位函數(shù)，將二階偏倒數(shù)運(yùn)算轉(zhuǎn)化為代數(shù)運(yùn)算，使求解過(guò)程更加簡(jiǎn)單，再運(yùn)用EXCEL作為計(jì)算工具得出結(jié)論的過(guò)程就更為簡(jiǎn)化。由于差分方程代替偏微分方程和差分方程本身的特點(diǎn)，在求解過(guò)程中存在一定的誤差，所以在計(jì)算過(guò)程中要根據(jù)具體問(wèn)題的精度要求選擇合適的網(wǎng)格劃分。而由于EXCEL里的單元格有限，并且在本文中用單元格描述真空罐的邊緣時(shí)用的是四舍五入取整的方法，這就使得其精度上受到限制，最終根據(jù)數(shù)據(jù)處理得到的數(shù)據(jù)繪制的電位圖不夠精細(xì)，但仍然可以為實(shí)際工程中的快速計(jì)算提供參考。

［1］ Giuliano V，Michele D S，F(xiàn)ranco G，Roberto B，Mario P.The INAF－IFSI Large Plasma Chamber［J］.INAF／IFSI，2010，（5）：1－5.

［2］ 謝處方.電磁場(chǎng)與電磁波（第4版）［M］.北京：高等教育出版社，2006.

［3］ 梁志輝.用差分法計(jì)算柱坐標(biāo)系的拉普拉斯方程［J］.信息系統(tǒng)工程，2009，（3）：104－105.

［4］ 盛劍霓.電磁場(chǎng)數(shù)值分析［M］.北京：科學(xué)出版社，1984.

［5］ 賈新民，嚴(yán) 文.有限差分法求解拉普拉斯方程［J］.昌吉學(xué)院學(xué)報(bào)，2009，（5）：105－109.

［6］ 趙宣銘，趙玉懷.典型數(shù)學(xué)物理方程的Excel數(shù)值解法［J］.紡織高?；A(chǔ)科學(xué)學(xué)報(bào)，2009，23（1）：118－126.

［7］ 何紅雨.電磁場(chǎng)數(shù)值計(jì)算法與 MATLAB實(shí)現(xiàn)［M］.武漢：華中科技大學(xué)出版社，2004.

［8］ 葛 超，王 蕾，曹秀爽.MATLAB技術(shù)大全［M］.北京：人民郵電出版社，2014.

［9］ 張雅男，徐 飛，葉 穎.MATLAB模擬靜電場(chǎng)與模擬靜電場(chǎng)實(shí)驗(yàn)的比較［J］.物理與工程，2008，（7）：32－33.