正六邊形電位分布的有限差分算法.doc

正六邊形電位分布的有限差分算法姓 名：指導老師： 學號： 目 錄正六邊形電位分布的有限差分算法5引言51有限差分法的基本原理步驟52正六邊形二維場域電位的有限差分算法62.1邊界處理62.2數(shù)學模型的建立72.3正六邊形電位分布的仿真程序93討論及分析13正六邊形電位分布的有限差分算法摘要：介紹了應(yīng)用有限差分法求電位分布的一般步驟，針對靜電場中軸對稱情形下的正六邊形場域的電位分布，建立了正三角形網(wǎng)格劃分的有限差分法的計算模型，給出了Matlab仿真的程序設(shè)計流程圖，并通過編程得到場域內(nèi)的電位分布圖形，對有限差分法的計算處理進行了討論和分析。關(guān)鍵詞：正六邊形；matlab；電位；有限差分法引言關(guān)于電磁場的數(shù)值計算，常用的方法有：有限差分法，時域有限差分法，有限單元法，矩量法，邊界元素法等。由于這些方法只能獲得近似解，因此，利用軟件進行仿真或者求數(shù)值解就顯得非常必要。本文基于Matlab，利用有限差分法求解靜電場中的正六邊形的電位分布。1有限差分法的基本原理步驟有限差分法的基本思想是把連續(xù)的定解區(qū)域用有限個離散點構(gòu)成的網(wǎng)格來代替，這些離散點稱作網(wǎng)格的節(jié)點；把連續(xù)的定解區(qū)域上的連續(xù)變量的函數(shù)用在網(wǎng)格上定義的離散變量函數(shù)來近似；把原方程和定解條件中的微商用差商來近似，積分用求和來近似，于是原微分方程和定解條件就近似地代之以代數(shù)方程組，即有限差分方程組，解此方程組就可以得到原問題在離散點上的近似解。然后利用差值方法便可以從離散解得到定解問題在整個區(qū)域上的近似解。有限差分法數(shù)值計算包括下列基本步驟：1)區(qū)域的離散或子區(qū)域劃分；2)插值函數(shù)的選擇；3)方程組的建立；4)求解方程組。2正六邊形二維場域電位的有限差分算法軸對稱電磁場問題是電工設(shè)備設(shè)計分析中常遇到的一大類問題。如圖1所示，邊長為b的正六邊形二維場域內(nèi)無電荷分布，6條邊上的電位（V）依次為1，-1，1，-1，1，-1，求場域內(nèi)的電位分布。2.1邊界處理由對稱性容易看出，正六邊形外接圓的D 條直徑EE，F(xiàn)F和GG 均為零電位線。因此，被這3條直徑切割成的6個正三角形區(qū)域的電位函數(shù)不獨立，而具有如下性質(zhì)。A點的電位與B點的電位滿足 (1)同理，有 (2) (3)其中，OA=OB=。當然，即使在-3030范圍內(nèi)，電位數(shù)據(jù)仍存在冗余現(xiàn)象。所以，本題的正六邊形二維場域電位分布的計算問題可以化為一個正三角形的電位分布問題，只要求出一個正三角形的電位分布，其它的就可以由（1）、（2）、（3）式的關(guān)系來確定。而一個正三角形場域中的電位計算，等價于下述拉普拉斯方程邊值問題。場域：=30和x=3b/2x03條直線圍成的等邊三角形區(qū)域，如圖1所示的OGE。邊界條件：=0（當=30）=1（當x= x0）但是，該問題是三角形場域，因此，如果用通常的正方形網(wǎng)格劃分邊界，那么邊界就不能恰好地落在網(wǎng)格上，這樣一方面給計算編程帶來麻煩，其次會使計算產(chǎn)生邊界取值的誤差。所以，針對場域形狀采用三角形網(wǎng)格劃分是處理邊界條件的好辦法。一般在進行網(wǎng)格劃分時采用對稱性網(wǎng)格形式，這樣既方便數(shù)學建模，也方便計算編程。 圖1 正六邊形場域的邊值問題 圖2 場域的正三角形網(wǎng)格劃分2.2數(shù)學模型的建立二維場域的拉普拉斯方程可以用有限差分法進行近似計算。首先把求解的區(qū)域劃分成網(wǎng)格，再把求解區(qū)域內(nèi)連續(xù)的場分布用求網(wǎng)絡(luò)節(jié)點上離散的數(shù)值解代替。網(wǎng)格必須劃分得充分細，才能達到足夠的精度。如圖2所示，對于正三角形場域OGE，采用正三角形網(wǎng)格劃分。其邊界全部由網(wǎng)格點來劃分，避免了邊界取值的誤差，也方便了計算編程。但域中任一點P的相鄰點有6個，因此，用有限差分法計算編程需另建數(shù)學模型。設(shè)每個正三角形網(wǎng)格邊長為a（稱為步長），網(wǎng)格節(jié)點（i，j）的電位為i,j，與其保持等距離的6個鄰點的電位分別為i,j+1，i,j-1，i-1,j，i-1,j-1，i+1,j，i+1,j+1。在a充分小的情況下，可以i,j為基點進行泰勒級數(shù)展開。 （4） （5） （6） （7） （8） （9）其中：為沿l方向的方向?qū)?shù)，為沿方向的方向?qū)?shù)。由于方向?qū)?shù)可表為所以，可得的二次方向?qū)?shù)為 （10）同理 （11）把（4）-（9）式相加，得把（10）、（11）式代入上式得 （12）對于（12）式，由于拉普拉斯方程為所以（12）式變?yōu)楸硎綼的4階無窮小，可以略去不計，則有限差分的數(shù)學表達式為 （13）2.3正六邊形電位分布的仿真程序?qū)D1網(wǎng)格節(jié)點數(shù)設(shè)置為3719=703，迭代精度為10-6。根據(jù)（13）式利用Matlab編制程序在計算機上運行，計算程序流程圖如圖3所示。程序運行得到計算結(jié)果的迭代次數(shù)為65。圖4為程序計算結(jié)果的圖示，它描述了正六邊形二維區(qū)域內(nèi)電位的等位線分布情況，其分布結(jié)果一目了然。圖3 計算程序流程圖1） 采用簡單迭代法求解簡單迭代法的特點是用前一次迭代得到的網(wǎng)絡(luò)點電位作為下一次迭代時的初值，迭代時計算公式為程序運行得到計算結(jié)果的迭代次數(shù)為65，最后電位數(shù)值解收斂于某一固定值。表1、表2列出了上、下兩個正三角形節(jié)點電位的差分運算結(jié)果。原則上，知道正六邊形二維場域里的任意一個正三角形中的電位分布，由（1）-（3）式可得到另外的5個正三角形的電位分布，但正六邊形二維場域里的6個正三角形的電位邊值畢竟有兩種：一種是外邊界值為1，一種是外邊界值為-1。 因此，只要把這兩種邊界的正三角形內(nèi)的電位值算出，剩下的4個正三角形內(nèi)電位分布就完全類似了，這樣就能直觀地分析正六邊形中的電位分布情況。本文以上、下正三角形為兩種不同邊值的場域為例，列出計算結(jié)果比較，在對角線上的電位值為0，相應(yīng)的對稱點的電位絕對值相等，符號相反，并沿著x軸越靠近中心，其電位絕對值越小，最后中心點O的電位也為零。表1 上三角形節(jié)點電位的差分運算結(jié)果0-1.00-1.00-1.00-1.00-1.00-1.00-1.00-1.00-1.00-1.00-1.00-1.00-1.00-1.00-1.00-1.00-1.0000-0.41-0.65-0.75-0.79-0.81-0.82-0.82-0.81-0.80-0.78-0.75-0.70-0.64-0.54-0.3900-0.28-0.46-0.57-0.62-0.65-0.65-0.64-0.62-0.58-0.53-0.46-0.37-0.2300-0.20-0.35-0.43-0.48-0.49-0.49-0.46-0.42-0.36-0.27-0.1600-0.15-0.26-0.32-0.35-0.35-0.32-0.28-0.21-0.1200-0.11-0.18-0.22-0.23-0.21-0.16-0.0900-0.08-0.21-0.13-0.11-0.0700-0.04-0.06-0.0400-1.0000表2 下三角形節(jié)點電位的差分運算結(jié)果001.00000.040.060.04000.080.210.130.110.07000.110.180.220.230.210.160.09000.150.260.320.350.350.320.280.210.12000.200.350.430.480.490.490.460.420.360.270.16000.280.460.570.620.650.650.640.620.580.530.460.370.23000.410.650.750.790.810.820.820.810.800.780.750.700.640.540.39001.001.001.001.001.001.001.001.001.001.001.001.001.001.001.001.001.000圖4為程序計算結(jié)果的圖示，它描述了正六邊形二維區(qū)域內(nèi)電位的等位線分布情況，其分布結(jié)果一目了然，與其物理分析的結(jié)果一致。源程序代碼見附錄。圖4 正六邊形二維區(qū)域內(nèi)的電位分布圖2） 采用超松弛法（SOR）求解簡單迭代法在解決問題時收斂速度比較慢，一般來說，實用價值不大。實際中常采用超松馳法， 相比之下它有兩點重大的改進。第一是計算每一網(wǎng)格點時，把剛才計算得到的鄰近點的電位新值代入，即在計算點（i，j）的電位時，把點（i-1，j）、點（i，j-1）和點（i-1，j-1）的電位用剛才算得的新值代入，即上式稱為松馳法或賽德爾法（relaxation method）。由于提前使用了新值，使得收斂速度加快。第二，再把上式寫成增量形式這時每次的增量（即上式右邊的第二項）就是要求方程局部達到平衡時應(yīng)補充的量。為了加快收斂，我們引進一個松馳因子w，將上式改寫為SOR法中w的取值對迭代公式的收斂速度影響很大，它的好壞直接影響到加速的快慢。為了保證迭代過程的收斂，必須要求0w2，超松弛法取1w2 。但是在1和2之間仍然有很多的取值，究竟如何取值沒有統(tǒng)一的規(guī)定，目前有學者提出了一些方法，如逐步實驗法、折半查找法、經(jīng)驗法和基于最小二乘法等。對于傳統(tǒng)的矩形網(wǎng)格劃分法，其最佳松弛因子為式中m、n為x、y方向的網(wǎng)格數(shù)。本文采用逐步實驗法找尋最優(yōu)松弛因子，首先采用如上的公式計算松弛因子得到w=1.7582，迭代次數(shù)k=88次，比采用普通迭代法時收斂速度要慢很多。可以得到此時w太大，逐步減小w的值，以得到最優(yōu)松弛因子。選取不同的w時，達到同樣的收斂精度10-6所需迭代次數(shù)如表3所示。表3 不同的w值所對應(yīng)的迭代次數(shù)w迭代次數(shù)w迭代次數(shù)1.7582881.6500571.5500391.4500291.4200271.3900271.3800281.3700291.3500311.300035由表3可以看到，當w為1.4左右時，收斂速度可以達到最快，迭代次數(shù)僅27次。3） 在該正六邊形場域內(nèi)加矩形導體為了與實際的電磁場問題結(jié)合起來，在該正六邊形場域內(nèi)加矩形導體，導體表面的電位（V）為1，計算此時的電位分布，其結(jié)果如圖5所示。圖5 正六邊形二維區(qū)域內(nèi)加矩形導體的電位分布圖3討論及分析（1）傳統(tǒng)的二維場域的電位數(shù)值差分計算，都是把整個場域網(wǎng)格劃分，再進行計算，而很少考慮場域的對稱性。而本文則是根據(jù)場域的對稱性，盡量把場域的計算區(qū)域變小，這樣可減少計算機的計算時間。對正六邊形場域處理就是把整個區(qū)域歸結(jié)為一個正三角形場域OGE的電位計算，其計算時間大大縮小。（2）傳統(tǒng)的二維場域網(wǎng)格劃分一般是矩形或正方形。本文采用了正三角形網(wǎng)格劃分，這是作者基于正六邊形場域所作的特殊處理，因為，這樣做可以使一些網(wǎng)格點落在邊界上，邊界網(wǎng)格點的值可準確確定，因此，可提高計算精度，由文計算結(jié)果可見，其正六邊形場域電位計算的精度是很高的。由此說明，為提高計算精度，除網(wǎng)格劃分變細之外，還可以采用特殊形狀網(wǎng)格劃分。（3）場域網(wǎng)格劃分的不同，其數(shù)值計算的數(shù)學模型也不同。傳統(tǒng)的矩形或正方形網(wǎng)格劃分便于數(shù)學建模及數(shù)值編程計算，本文采用正三角形網(wǎng)格劃分，對數(shù)學建模及數(shù)值編程帶來一定的麻煩，因此，在方向?qū)?shù)概念的基礎(chǔ)上建立了某點相鄰電位之和的1/6模型，這是與傳統(tǒng)的四方形網(wǎng)格數(shù)值計算模型所不同的，并且其截斷誤差在4階無窮小層次上，具有比較高的精度。（4）迭代方法的不同，會帶來收斂速度的不同。本文采用了超松弛迭代法，討論了不同的松弛因子的選擇對收斂速度的影響，并與普通的迭代方法進行了比