电磁场数值计算上机作业

电磁场数值计算上机作业报告一、有限差分法及原理有限差分法基本思想是把连续的定解区域用有限个离散点构成的网格来代替,这些离散点称作网格的节点；把连续定解区域上的连续变量的函数用在网格上定义的离散变量函数来近似；把原方程和定解条件中的微商用差商来近似,积分用积分和来近似,于是原微分方程和定解条件就近似地代之以代数方程组,即有限差分方程组,解此方程组就可以得到原问题在离散点上的近似解.然后再利用插值方法便可以从离散解得到定解问题在整个区域上的近似解.在采用数值计算方法求解偏微分方程时,若将每一处导数由有限差分近似公式替代,从而把求解偏微分方程的问题转换成求解代数方程的问题。求解拉普拉斯方程：22220xy为简单起见，将场域分成足够小的正方形网格，网格线之间的距离为h，0h。节点0、1、2、3、4上的电位分别用0、1、2、3和4表示。点1、点3在x0处可微，沿x方向在x0处的泰勒级数展开式为2323100002311()()()()().2!3!hhhxxx2323300002311()()()()().2!3!hhhxxx点2、点4在y0处可微，沿y方向在y0处的泰勒级数展开式为2323200002311()()().2!3!hhhyyy2323400002311()()()()().2!3!hhhyyy忽略高次项22212340000224()()()4hxy稍作变化得到拉普拉斯方程的五点差分格式：123404利用超松弛迭代法求解以上差分方程，二维场拉普拉斯方程等距剖分差分格式公式为：𝜑ⅈ,𝑗𝑛+1=𝜑ⅈ,𝑗𝑛+𝛼4(𝜑𝑖+1,𝑗𝑛+𝜑ⅈ,𝑗+1𝑛+𝜑ⅈ−1,𝑗𝑛+1+𝜑ⅈ,𝑗−1𝑛+1−4𝜑ⅈ,𝑗𝑛)其中𝛼为超松驰因子，1≤𝛼＜2.计算流程如图1所示：图1超松驰迭代法计算流程利用上述方法求解两题。二、第一题计算长直接地金属槽中的电场分布。金属槽横截面如图所示，其侧壁与底面电位均为零，顶盖电位相对值为10。槽内电位函数满足拉普拉斯方程。计算槽内电位分布。要求：（1）先用正方形网格粗分，每边取4个网格计算，取不同的松弛因子，比较其收敛速度。取计算精度为千分之一。（2）划分网格加倍，计算电位分布，并与上面计算结果比较。（1）用正方形网格，每边取4个网格计算。图2题1示意图网格分割如图3所示：取α=1，计算得结果如图4：图3每边4个网格分割图44个网格分割计算结果具体数值表示结果如图5所示：取不同的超松驰因子α，得不同的收敛速度如下表：超松驰因子α11.11.21.31.41.51.61.71.81.9循环次数1411810121619274482可见随着α从1到2的增加，收敛速度先变快后变慢。在α=1.2时收敛速度最快。图54个网格分割计算结果数值表示表1不同超松驰因子α下的迭代次数(2)划分网格加倍，计算结果如下图所示图68个网格分割计算结果具体数值表示为：与划分为4个网格的结果做比较有下表：X轴Y轴0.000.25a0.5a0.75aa8格4格8格4格8格4格8格4格8格4格a10.0010.0010.0010.0010.0010.0010.0010.0010.0010.000.75a0.000.004.314.295.365.274.314.290.000.000.5a0.000.001.841.872.502.501.841.870.000.000.25a0.000.000.690.710.960.980.690.710.000.000.000.000.000.000.000.000.000.000.000.000.00可见在拟合过程中，划分网格较少时，距离y=a较近处，数据偏小，而距离y=0较近处，数据偏大。可见划分网格较多时可以使得到的电势分布趋于平滑。表2相同节点位置不同网格划分计算数值对比图78个网格分割计算结果数值表示三、第二题计算电机电枢槽气隙磁位分布：如图所示，忽略曲度效应，采用直角坐标系，计算二维恒定场。令定子和转子磁铁中μ=∞，EF是对称线，考虑对称性有：𝜙1|𝐴𝐹=100𝜙2|𝐵𝐶𝐷𝐸=0|𝐴𝐵=0|𝐸𝐹=0求G1（ABCDEF）内部和边界上全部节点磁位，计算精度取ω≤1/1000。设g=8，h=66，BC=6，DE=5。已知磁位𝜙也满足拉普拉斯方程，在所求区域内部，利用超松弛迭代法给出二维场拉普拉斯方程等距剖分差分格式公式为：𝜑ⅈ,𝑗𝑛+1=𝜑ⅈ,𝑗𝑛+𝛼4(𝜑𝑖+1,𝑗𝑛+𝜑ⅈ,𝑗+1𝑛+𝜑ⅈ−1,𝑗𝑛+1+𝜑ⅈ,𝑗−1𝑛+1−4𝜑ⅈ,𝑗𝑛)在AB处，有|𝐴𝐵=0，转化为差分形式可得𝜑𝑖+1,𝑗𝑛=𝜑ⅈ−1,𝑗𝑛+1，所以在AB边界处的差分形式方程为：𝜑ⅈ,𝑗𝑛+1=𝜑ⅈ,𝑗𝑛+𝛼4(2×𝜑𝑖+1,𝑗𝑛+𝜑ⅈ,𝑗+1𝑛+𝜑ⅈ,𝑗−1𝑛+1−4𝜑ⅈ,𝑗𝑛)同理在EF处满足的差分形式方程为：图8题2示意图𝜑ⅈ,𝑗𝑛+1=𝜑ⅈ,𝑗𝑛+𝛼4(2×𝜑𝑖−1,𝑗𝑛+𝜑ⅈ,𝑗+1𝑛+𝜑ⅈ,𝑗−1𝑛+1−4𝜑ⅈ,𝑗𝑛)按正方形网格划分，网格大小为1×1，根据对称性，得整个区域的磁位如图9：在G1区域中，在每个节点用数值表示为（放大可看具体数值）：图9题2计算结果至此求得全部磁位。附：matlab代码第一题：clc;x=5;y=5;foralpha_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;fori=2:(x-1)forj=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&&k10000)Fai_Last=Fai_Now;fori=2:(x-1)forj=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=k+1;endk_test(alpha_test-9)=k;endlx=linspace(0,a,x);ly=linspace(0,a,y);figure(1);surf(lx,ly,Fai_Now);xlabel('x轴/a','FontName','微软雅黑','FontSize',20);ylabel('y轴/a','FontName','微软雅黑','FontSize',20);zlabel('z轴/V','FontName','微软雅黑','FontSize',20);图10题2计算结果数值表示title('电势φ分布','FontName','微软雅黑','FontSize',20);set(gca,'XTick',0:1/(x-1):1,'FontSize',12);set(gca,'YTick',0:1/(y-1):1,'FontSize',12);holdonfigure(2);[x_s,y_s]=meshgrid(0:x-1);x_s=x_s./(x-1);y_s=y_s./(y-1);plot(x_s,y_s,'k',x_s',y_s','k');set(gca,'XTick',0:1/(x-1):1,'FontSize',12);set(gca,'YTick',0:1/(y-1):1,'FontSize',12);xlabel('x轴/a','FontName','微软雅黑','FontSize',20);ylabel('y轴/a','FontName','微软雅黑','FontSize',20);title('电势φ分布','FontName','微软雅黑','FontSize',20);holdonfori=1:xforj=1:ytext_show=strcat('φ=',num2str(Fai_Now(i,j)));text((i-1)/(x-1),(j-1)/(y-1),text_show,'color','r','FontSize',13);endendholdoffxlswrite('homework.xls',rot90(Fai_Now));第二题：clc;h=66;g=8;BC=6;DE=5;length_x=BC+DE;length_y=h+g;div_x=length_x+1;div_y=length_y+1;alpha_test=10;Fai_Last=zeros(div_x,div_y);Fai_Now=zeros(div_x,div_y);Fai_Last(:,div_y)=100;Fai_Now(:,div_y)=100;a=1;alpha=alpha_test/10;fori=1:(div_x)forj=1:(div_y)if(i1&&idiv_x&&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));endif(iBC+1&&idiv_x&&j1&&j=h+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));endif(i==div_x&&j1&&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));endif(i==1&&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));endendendk=1;while(max(max(abs((Fai_Last-Fai_Now)./Fai_Last)))0.001&&k10000)Fai_Last=Fai_Now;fori=1:(div_x)forj=1:(div_y)if(i1&&idiv_x&&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));endif(iBC+1&&idiv_x&&j1&&j=h+1)Fai_Now(i,j)=Fai_Last(i,j)+alpha/4*(Fai_Last(i,j+1)+