恒流导体表面电荷分布形成的计算机模拟

恒流导体表面电荷分布形成的计算机模拟

项长铖
中国科学技术大学09级物理2班Email:willon@mail.ustc.edu.cn

摘要


     本文将建立二维模型通过计算机模拟的方法进行恒流导体表面电荷分布的分析,并对模型的精确性进行了初步讨论。

引入

     在任何一本电磁学教科书上都会介绍稳恒电流的稳恒条件为\oint_{D} \vec{j}\cdot d\vec{S}=0\nabla\cdot\vec{j}=0 。根据这一结论,我们将得到一个十分有趣的结论:如左图所示,在导体和外界的分界面上取高度\Delta{h}的柱面D为高斯面,当\Delta{h}趋于0时应用上述稳恒条件\oint_{D} \vec{j}\cdot d\vec{S}\approx (\vec{j}\cdot\vec{n})S=0,可以看到\vec{j}始终与导体界面的法向量\vec{n}垂直,也就是说导体表面正是电流管的表面,又根据欧姆定律\vec{j}=\sigma\vec{E}导体表面又是电场线管的表面。更有意思的是,无论导体被弯折成何种形状,电场线总是被“束缚”在导体里。然而,刚通电时空间的电场分布如右图所示,是呈发散状的,这与稳恒后的情形不同。
是什么将电场“束缚”在导体中呢?由于导体内部无电荷分布,自然是导体表面的电荷了。那么这种电荷究竟如何分布呢?传统方法是先求解混合边界条件(即导体两端电势固定,其他表面则满足Neuman条件 \frac{\par{\varphi}}{\par{n}} =0)的Laplace方程\nab^{2}{\varphi}=0得到电势分布\varphi,而后\vec{E}=-\nab{\varphi}得出外部电场分布,最后由\sig=\vec{E}\cdot\vec{n}得到电荷分布。
     这种方法具有很大的局限性,一旦边界形状稍有复杂,就很难求解Laplace方程。但是得益于计算机数值计算技术的发展,我们可以通过计算机模拟表面电荷分布形成的过程并得到最终的电荷分布。

 


第一部分 模型的建立

为了讨论的简便,这里只研究二维的模型,三维模型完全可以用一样的方法处理。事实上,二维模型也对应了三维模型中的某些情况。

     计算机数值模拟需要的是离散的模型。建立模型的第一步是将实际的电磁场连续域离散化,即取有限空间中离散的网格点阵(设相邻格点的间距为很小的h),用这些点的离散的电势场值\varphi(i,j)和电荷量q(i,j)逼近连续域中对应的某一极小体积空间处的真实值,同时用差商\frac{f(x+h,y)-f(x,y)}{h}代替微商\frac{\par{f}}{\par{x}}=\lim_{h \rightarrow 
0}{\frac{f(x+h,y)-f(x,y)}{h}} ,最终将泊松方程\frac{\par^2\varphi}{\par{x^2}}+\frac{\par^2\varphi}{\par{y^2}}=\rho变为差分形式\varphi(i+1,j)+\varphi(i-1,j)+\varphi(i,j+1)+\varphi(i,j-1)-4\varphi(i,j)=q(i,j)。根据给出的边界条件:空间边界电势为0(若空间取得足够大可以近似认为该空间的边界为无穷远因而电势为0。当然即使这个空间不一定很大,仍然可以假想是吧导体放入了一个与给定空间边界重合的接地导体壳中,同样由边界电势为0)再加上差分形式的泊松方程\varphi(i+1,j)+\varphi(i-1,j)+\varphi(i,j+1)+\varphi(i,j-1)-4\varphi(i,j)=q(i,j)^{[1]} (1)就可以解出某一时刻的电势分布了。然后再用差商替代微商的方法求出E_{x}(i,j)=-\frac{\varphi(i,j+1)-\varphi(i,j)}{h},E_{y}(i,j)=-\frac{\varphi(i+1,j)-\varphi(i,j)}{h}(2)得出电场分布,最后是\vec{j}=\sig\vec{E}(方便计算取\sig=1)求出电流分布。

     第一步解决的仅是将某一是时间片断中电场、电流场域的离散化,而要处理的问题是时变的,故第二步是将导体内电流的行为离散化。实际的过程是当\Del{t},\Del{h}\rightarrow {0}时,\Del{q(x,y)}=(\nab\cdot{\vec{j}(x,y)})\Del{t}\Del{h}^{2},还是利用差商替代微商的方法,\nab\cdot{\vec{j}}=\frac{\par{j_{x}}}{\par{x}}+\frac{\par{j_{y}}}{\par{y}}=\sig(\frac{\par{E_{x}}}{\par{x}}+\frac{\par{E_{y}}}{\par{y}})=\frac{[(E_{x}(i,j+1)-E_{x}(i,j)+(E_{y}(i+1,j)-E_{y}(i,j))]}{\Del{h}}(若有某一格不是导体则电场分量用0代替)则\Del{q(i,j)}=[(E_{x}(i,j+1)-E_{x}(i,j)+(E_{y}(i+1,j)-E_{y}(i,j))]\Del{h}\Del{t}(3)(同样为计算方便取\Del{h}=1)(以下是一些实际问题:在真实情况下导体内电荷变化应为0,但是由以上公式,只能证明\Del{q}=O(\Del{h}^2)O(\Del{t})虽然很小但是这将导致\frac{d\rho}{dt}=\nab\cdot\vec{j}=\sig\nab\cdot\vec{E}=\rho\ne 0从而使电荷逐渐积累,因此在处理时只考虑导体边界上的电荷变化(**)

根据这些我们就能画出流程图并通过matlab编程来求解了。流程图如下:



第二部分 对模拟结果的讨论

     根据图2直导体刚通电时径向电场从两端到中间逐渐减小,粗略地看,要将电场“束缚”在导体中,表面电荷分步应该是两端多于中间。以下是模拟结果:



模拟结果与确实预想一致。

再考虑弯折的导线,要将电场随导线“弯折”,不可能再向上面一样中间电荷分布极少,很可能出现另一峰值。模拟结果如下:





和预料的相差无几,还得到另一结果:折角内部的电荷分布比外部相对较多,且并未明显呈现中段迅速减少的现象,出于将电场“束缚”在导体内的最终目的,这也是可以理解的。
仔细观察模拟结果仍有不够令人满意的地方:在电势分布图中,导体的负极区附近出现了一个“空腔”,也就是说电荷分布有一个不正常的变化趋势导致该处电场为0。猜测是因(**)处考虑的原因,由于\Del{t}不够小导致一开始的电荷变化不满足电荷守恒,从而造成这一现象。修改程序中的\Del{t}后,确有改观,但限于个人计算机计算的水平有限,仍不能将其修改到足够小以使“空腔”消失,需要寻找运算能力更高的计算机,或者寻找更合理的模型才能解决这一缺陷。


  本文结论

     该模型能够初步验证关于导体表面分布的一些猜想,同时也能得出一些较可靠的结论:如,直导体表面电荷分步应该是两端多于中间;弯折导体还满足弯折角内部电荷多于外部,并在弯折处电荷分布有一峰值。但是,由于模型一些地方存在缺陷,要完全靠该模型模拟实际情况还需要有强大计算能力的计算机,因此该模型有待改进,或者由更精确的模型代替它。

      参考文献:[1]何红雨著,电磁场数值计算法与MATLAB实现,华中科技大学出版社,2005
                  [2]胡友秋 程福臻 叶邦角,电磁学与电动力学(上) ,科学出版社,200


Follow @Micaiah_Yin