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

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