电磁场数值计算上机作业.doc

.电磁场数值计算上机作业报告一、 有限差分法及原理有限差分法基本思想是把连续的定解区域用有限个离散点构成的网格来代替,这些离散点称作网格的节点；把连续定解区域上的连续变量的函数用在网格上定义的离散变量函数来近似；把原方程和定解条件中的微商用差商来近似,积分用积分和来近似,于是原微分方程和定解条件就近似地代之以代数方程组,即有限差分方程组 ,解此方程组就可以得到原问题在离散点上的近似解.然后再利用插值方法便可以从离散解得到定解问题在整个区域上的近似解.在采用数值计算方法求解偏微分方程时,若将每一处导数由有限差分近似公式替代,从而把求解偏微分方程的问题转换成求解代数方程的问题。求解拉普拉斯方程：为简单起见，将场域分成足够小的正方形网格，网格线之间的距离为h，。节点0、1、2、3、4上的电位分别用、和表示。点1、点3在x0处可微，沿x方向在x0处的泰勒级数展开式为点2、点4在y0处可微，沿y方向在y0处的泰勒级数展开式为忽略高次项稍作变化得到拉普拉斯方程的五点差分格式：利用超松弛迭代法求解以上差分方程，二维场拉普拉斯方程等距剖分差分格式公式为：,jn+1=,jn+4i+1,jn+,j+1n+-1,jn+1+,j-1n+1-4,jn其中为超松驰因子，12.计算流程如图1所示：图1 超松驰迭代法计算流程利用上述方法求解两题。二、 第一题计算长直接地金属槽中的电场分布。金属槽横截面如图所示，其侧壁与底面电位均为零，顶盖电位相对值为10。槽内电位函数满足拉普拉斯方程。计算槽内电位分布。图2 题1示意图要求：（1）先用正方形网格粗分，每边取4个网格计算，取不同的松弛因子，比较其收敛速度。取计算精度为千分之一。（2）划分网格加倍，计算电位分布，并与上面计算结果比较。（1） 用正方形网格，每边取4个网格计算。图3 每边4个网格分割网格分割如图3所示：图4 4个网格分割计算结果取=1，计算得结果如图4：图5 4个网格分割计算结果数值表示具体数值表示结果如图5所示：表1 不同超松驰因子下的迭代次数取不同的超松驰因子，得不同的收敛速度如下表：超松驰因子11.11.21.31.41.51.61.71.81.9循环次数1411810121619274482可见随着从1到2的增加，收敛速度先变快后变慢。在=1.2时收敛速度最快。图6 8个网格分割计算结果(2)划分网格加倍，计算结果如下图所示图7 8个网格分割计算结果数值表示具体数值表示为：与划分为4个网格的结果做比较有下表：X轴Y轴0.00 0.25a0.5a0.75aa8格4格8格4格8格4格8格4格8格4格a10.00 10.00 10.00 10.00 10.00 10.00 10.00 10.00 10.00 10.00 0.75a0.00 0.00 4.31 4.29 5.36 5.27 4.31 4.29 0.00 0.00 0.5a0.00 0.00 1.84 1.87 2.50 2.50 1.84 1.87 0.00 0.00 0.25a0.00 0.00 0.69 0.71 0.96 0.98 0.69 0.71 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 表2 相同节点位置不同网格划分计算数值对比可见在拟合过程中，划分网格较少时，距离y=a较近处，数据偏小，而距离y=0较近处，数据偏大。可见划分网格较多时可以使得到的电势分布趋于平滑。三、 第二题图8 题2示意图计算电机电枢槽气隙磁位分布：如图所示，忽略曲度效应，采用直角坐标系，计算二维恒定场。令定子和转子磁铁中=，EF是对称线，考虑对称性有：1AF=100 2BCDE=0xAB=0 xEF=0求G1（ABCDEF）内部和边界上全部节点磁位，计算精度取1/1000。设g=8，h=66，BC=6，DE=5。已知磁位也满足拉普拉斯方程，在所求区域内部，利用超松弛迭代法给出二维场拉普拉斯方程等距剖分差分格式公式为：,jn+1=,jn+4i+1,jn+,j+1n+-1,jn+1+,j-1n+1-4,jn在AB处，有 xAB=0 ，转化为差分形式可得i+1,jn=-1,jn+1 ，所以在AB边界处的差分形式方程为：,jn+1=,jn+42i+1,jn+,j+1n+,j-1n+1-4,jn同理在EF处满足的差分形式方程为：,jn+1=,jn+42i-1,jn+,j+1n+,j-1n+1-4,jn图9 题2计算结果按正方形网格划分，网格大小为11，根据对称性，得整个区域的磁位如图9：在G1区域中，在每个节点用数值表示为（放大可看具体数值）：图10 题2计算结果数值表示至此求得全部磁位。附：matlab代码第一题：clc;x = 5; y = 5;for alpha_test=10:19;Fai_Last = zeros(x,y);Fai_Now = zeros(x,y);Fai_Last(:,y) = 10;Fai_Now(:,y) = 10;a = 1;alpha = alpha_test/10;for i=2:(x-1) for j=2:(y-1) Fai_Now(i,j) = Fai_Last(i,j)+alpha/4*(Fai_Last(i,j+1)+Fai_Last(i+1,j). +Fai_Now(i,j-1)+Fai_Now(i-1,j)-4*Fai_Last(i,j); endendk=1;while(max(max(abs(Fai_Last-Fai_Now)./ Fai_Last)0.001&k1&ih+1&jBC+1&i1&j1&jh+1&j0.001&k1&ih+1&jBC+1&i1&j1&jh+1&jdiv_y) Fai_Now(i,j) = Fai_Last(i,j)+alpha/4*(Fai_Last(i,j+1)+Fai_Last(i+1,j). +Fai_Now(i,j-1)+Fai_Now(i+1,j)-4*Fai_Last(i,j); end endend k=k+1;endk_test(alpha_test-9)=k;Fai_Final = zeros(BC+DE)*2+1,h+g+1);Fai_Final(1:BC+DE+1,1:h+g+1) = Fai_Now;Fai_temp = fliplr(Fai_Now);Fai_Final(BC+DE+1):(BC+DE)*2+1,1:h+g+1)=Fai_temp;lx = linspace(0,(BC+DE)*2+1,(BC+DE)*2+1);ly = linspace(0,h+g+1,h+g+1);figure(1);colormap(jet);surf(lx,ly,Fai_Final);set(gca,XTick,1:(BC+DE)*2+1,FontSize,3);set(gca,YTick,1:5:h+g+1,FontSize,3);xlabel(x轴 ,FontName,微软雅黑,FontSize,20);ylabel(y轴 ,FontName,微软雅黑,FontSize,20);zlabel(z轴 ,FontName,微软雅黑,FontSize,20);title(磁标势分布,FontName,微软雅黑,FontSize,20);hold on figure(2);x_s,y_s=meshgrid(0:(BC+DE)*2+1,0:h+g+1);plot(x_s,y_s,k,x_s,y_s,k);xlabel(x轴,FontName,微软雅黑,FontSize,20);ylabel(y轴,FontName,微软雅黑,FontSize,20);title(磁标势分布,FontName,微软雅黑,FontSize,20);axis(0 (B