靜電場(chǎng)電荷面密度與曲率關(guān)系的Matlab模擬研究

陳曙英，陳 鋼，崔宏基

(浙江工業(yè)大學(xué)應(yīng)用物理系，浙江 杭州 310023)

1 引言

對(duì)于避雷針并不陌生，它主要是利用尖端放電的原理引導(dǎo)雷電放電，再通過(guò)接地引線(xiàn)和裝置將雷電引入大地，從而達(dá)到保護(hù)建筑的目的。通過(guò)對(duì)尖端放電現(xiàn)象的研究，人們發(fā)現(xiàn)隨著導(dǎo)體表面曲率的增大，該位置的電荷密度也隨之升高。在強(qiáng)電場(chǎng)的作用下，導(dǎo)體中比較尖銳的部分則會(huì)發(fā)生放電現(xiàn)象，這就是尖端放電。

尖端放電原理的應(yīng)用有很多，如工業(yè)煙囪的除塵裝置、電子打火裝置等[1]。雖然尖端放電現(xiàn)象的應(yīng)用很廣泛，但是大多數(shù)電磁類(lèi)教科書(shū)中卻很少具體的討論表面電荷密度與導(dǎo)體曲率之間的關(guān)系。通過(guò)相關(guān)的數(shù)值模擬，分析并得到電荷面密度與曲率的關(guān)系，可以更好的豐富這些方面的研究。

2 理論分析

在講電荷面密度與曲率的關(guān)系前，必須先了解主曲率和平均曲率的概念。在已知曲面上取一點(diǎn)E，并且曲面在E點(diǎn)的法線(xiàn)為z軸，通過(guò)z軸可以得到許多剖切平面，如圖1所示，向上為z軸。曲面和每個(gè)剖切平面都相交，他們的交線(xiàn)是一條平面曲線(xiàn)，每條平面曲線(xiàn)在E點(diǎn)都有一個(gè)曲率半徑。不一樣的平面曲線(xiàn)在E點(diǎn)的曲率半徑一般是不相等的。這些曲率半徑中，找出最小和最大的曲率半徑，把這兩個(gè)稱(chēng)之為主曲率半徑，分別用R1和R2表示[2]。

圖1 曲率示意圖

在眾多的曲面論的問(wèn)題中，起決定性作用一般都不是主曲率本身，而是由它們的平均值和乘積來(lái)決定的，即平均曲率和高斯曲率。對(duì)于一些簡(jiǎn)單的曲面，可以直接通過(guò)定義求出平均曲率和高斯曲率[3]。如式（1）所示：

H就稱(chēng)為曲面在給定點(diǎn)E 的平均曲率，K就是E 點(diǎn)的總曲率或高斯曲率。

本文主要是利用MATLAB 模擬計(jì)算出相關(guān)曲面的曲率和其表面的電荷面密度，然后利用數(shù)值計(jì)算找出曲率與電荷面密度之間的關(guān)系。對(duì)于曲面，可以選擇橢球面、拋物面和雙曲面。并且經(jīng)過(guò)相關(guān)的公式推導(dǎo)，可以得到三個(gè)曲面的曲率和電荷面密度的嚴(yán)格數(shù)學(xué)關(guān)系，其中式（2）為橢球面，式（3）為拋物面，式（4）為雙曲面。式（2）中K 為總曲率，Q 為橢球?qū)w的總電荷，a和b代表赤道半徑，c為極半徑；式（3）中σ0為拋物面頂點(diǎn)（，0，0）處電荷面密度，k為常數(shù)；式（4）中a為實(shí)長(zhǎng)軸，b為實(shí)短軸，σ0為雙曲面頂點(diǎn)（，0，0）上的電荷面密度，k為常數(shù)。

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

計(jì)算中選擇的導(dǎo)體的曲面為橢球面、拋物面和雙曲面，并且給定整個(gè)導(dǎo)體所帶電荷數(shù)，然后利用MATLAB中的pdetool 工具模擬出電場(chǎng)的空間分布圖。其中橢球?qū)w的電場(chǎng)分布如圖2 所示。為了得到電荷面密度與曲率的關(guān)系，必須取得導(dǎo)體表面的電荷面密度，但pdetool 工具并不能直接給出電荷面密度的值。由于真空中靜電場(chǎng)中導(dǎo)體表面附近的場(chǎng)強(qiáng)滿(mǎn)足：

而導(dǎo)體表面附近任一點(diǎn)的E 值已經(jīng)通過(guò)數(shù)值模擬計(jì)算得到，因此可以直接換算出該點(diǎn)的電荷面密度的值，即，可以在導(dǎo)體表面附近的空間取得若干個(gè)點(diǎn)，然后依次得到各個(gè)點(diǎn)的電荷面密度。

對(duì)于橢球面導(dǎo)體，由于橢圓的長(zhǎng)短軸具有特殊性，因此選擇點(diǎn)1（0，0.3）和點(diǎn)2（0.5，0）兩個(gè)點(diǎn)。從圖2 中得到E1=5.31，E2=6.75。再根據(jù)σ=Eε0，可以得出兩點(diǎn)的電荷面密度σ1和σ2。同時(shí)，利用數(shù)學(xué)工具，可以得到兩點(diǎn)的K值，分別為K1=1.20，K2=5.56。根據(jù)公式（5），E與σ成正比，所以，而。

圖2 橢球面導(dǎo)體電場(chǎng)分布圖

對(duì)于橢球面電荷密度與曲率四分之一次方成正比[5]，而數(shù)值模擬給出的比值約為0.68，雖然比值之間存在誤差，但還是在誤差范圍之內(nèi)。對(duì)于造成誤差的原因主要有以下幾種：首先是數(shù)值模擬的橢圓不是真正意義上的橢圓，具有棱角；其次有限元計(jì)算時(shí)空間要進(jìn)行是三角分割，對(duì)計(jì)算有一定影響。

對(duì)于拋物線(xiàn)，在取點(diǎn)之前需要先確定定義域，然后任取兩個(gè)坐標(biāo)點(diǎn)。分別以?xún)蓚€(gè)坐標(biāo)點(diǎn)為中心，得出附近點(diǎn)的具體坐標(biāo)值和對(duì)應(yīng)的電壓值。通過(guò)數(shù)據(jù)的篩選，選擇拋物線(xiàn)的頂點(diǎn)坐標(biāo)（0，0），并且隨機(jī)取得坐標(biāo)（0.3，0.09），得出兩個(gè)點(diǎn)曲率之比為，數(shù)值計(jì)算得到的兩點(diǎn)的電荷面曲率之比為，數(shù)值計(jì)算和理論分析符合的非常好。

對(duì)于雙曲線(xiàn)，選擇的模型為（y+1)2-x=1，定義域?yàn)閇-1，1]，得到雙曲線(xiàn)模擬坐標(biāo)三角形分割圖，如圖2。對(duì)于坐標(biāo)點(diǎn)的選擇，選擇雙曲線(xiàn)下半部分的頂點(diǎn)（0，0）和上半部分的一點(diǎn)（0.30，0.04）。經(jīng)過(guò)計(jì)算，得出兩點(diǎn)電荷面密度之比為，理論分析給出的曲率的四分之一次方之比為。通過(guò)不同坐標(biāo)點(diǎn)的電荷面密度和曲率的比值，發(fā)現(xiàn)公式（5）仍然得到滿(mǎn)足。通過(guò)橢圓模型、雙曲線(xiàn)模型和拋物線(xiàn)模型的建立，證明了在誤差允許的范圍內(nèi)，橢圓、拋物線(xiàn)和雙曲線(xiàn)表面的電荷密度均與該處曲率的四分之一成正比。

圖3

之前的實(shí)驗(yàn)中，都是選取若干個(gè)點(diǎn)進(jìn)行比較，雖然最后得出了電荷面密度與曲率的關(guān)系，但說(shuō)服力不足，因此選取整條邊來(lái)進(jìn)行證明。通過(guò)數(shù)據(jù)的處理，得出拋物線(xiàn)、雙曲線(xiàn)和橢圓的一整條邊上的E值和曲率。通過(guò)Matlab 中的Cftool 工具，擬合出E值與曲率的關(guān)系圖如下：

圖4 橢圓擬合曲線(xiàn)

圖5 雙曲線(xiàn)擬合曲線(xiàn)

圖6 拋物線(xiàn)擬合曲線(xiàn)

通過(guò)擬合，得出擬合函數(shù)為E=a×Kb。其中，關(guān)于a和b 的取值，橢圓為a=3.847（3.772，3.923），b=0.25；雙曲線(xiàn)為a=3.584（3.373，3.795），b=0.25；拋物線(xiàn)為a=3.754（3.659，3.85），b=0.25。通過(guò)擬合曲線(xiàn)可以看出，E和K關(guān)系基本滿(mǎn)足E與K的四分之一次方成正比，所以可以推導(dǎo)出二次曲面導(dǎo)體表面電荷密度與曲率的四分之一次方成正比這一結(jié)論。

4 結(jié)論

利用軟件模擬出靜電感應(yīng)過(guò)程，并從中得出一系列坐標(biāo)點(diǎn)以及該點(diǎn)的電壓值。經(jīng)過(guò)數(shù)據(jù)的處理和計(jì)算，可以得到對(duì)應(yīng)點(diǎn)的電荷面密度及其表面的對(duì)應(yīng)曲率。將任意兩點(diǎn)之間的數(shù)據(jù)進(jìn)行對(duì)比，可以得出電荷面密度與曲率的對(duì)應(yīng)關(guān)系。通過(guò)本實(shí)驗(yàn)的結(jié)果，證明了通過(guò)軟件的模擬是可以更加方便的證明，對(duì)于二次曲面，電荷面密度與曲率的四分之一次方成正比。