电磁场数值计算作业报告

《电磁场数值计算作业报告》由会员上传分享，免费在线阅读，更多相关内容在行业资料-天天文库。

1、《电磁场数值计算》—有限元法报告第一题（一）问题描述及数学物理模型的建立①①9④5⑤6（12）128③⑥（11）11②3⑦4⑩7①1⑧2⑨10有一矩形场区，尺寸为6×9，如图1所示，在区域内部J=0，的左边界A=0,右边界A=100，上下边界满足：，媒质均为的磁导率为，利用有限元法求A的分布首先，此问题的微分方程及其定解条件为：图1求解场域及网格划分如果利用有限元来求解，那么对应的泛函为（注：分别为上边界和下边界）：第一类边界条件为：（二）场域的离散化及单元节点的编号处理网格的划分如图1所示，先求出每个单元的刚度矩阵，然后再进行整体合成，得到系数矩阵，具体的做法可参考教材，

2、注意，这里左边界和右边界属于第一类边界，因此，利用下式进行计算：这里，因为J=0,所以=0，节点编号如图1所示，这里，由于网格的划分和节点的编号都不是很规则，但节点数目和网格数目不是很多，所以，在程序中，通过穷举的办法，找到相应单元的（i，j，m）的值，以及相应的节点坐标。各单元的（i,j,m）（按逆时针方向）如下：1(3,7,1)2(3,8,7)3(9,8,3)4(5,9,3)5(6,5,3)6(6,3,4)7(4,3,1)8(4,1,2)9(4,2,10)10(11,4,10)11(12,4,11)12(12,6,4)节点i的坐标为：ifi=1,3,5thenx(i)=

3、3,y(i)=1.5(i-1)ifi=2,4,6thenx(i)=6,y(i)=1.5iifi=7,8,9thenx(i)=0,y(i)=3(i-7)ifi=10,11,12thenx(i)=9,y(i)=3(i-10)边界节点及其位函数的值为：7（0）8（0）9（0）10（100）11（100）12（100）然后建立有限元方程组，进行求解，方程组的建立过程参考教材（三）利用数值模拟计算按照有限元法计算的相关步骤，用FORTRAN90语言编写的程序清单如下：programhomework3_6implicitnoneintegeri,j,m,k,e,countreal*8s

4、,new,oldreal*8winteger::z(3)real*8::x(12),y(12),b(12),c(12)real*8::a(12,12),r(6),d(7:12),u(6)doi=7,9!x(i),y(i)为节点坐标x(i)=0y(i)=3*(i-7)d(i)=0!d(i)为第一类边界上的节点的磁位值enddodoi=10,12x(i)=9y(i)=3*(i-10)d(i)=100enddodoi=1,5,2x(i)=3y(i)=3*(i-1)/2enddodoi=2,6,2x(i)=6y(i)=3*(i-2)/2enddodoi=1,12doj=1,12a(

5、i,j)=0!系数矩阵的初始化enddoenddos=0.5*3*3dok=1,12!k为单元的标号callbianhao(k,z)!调用子例行程序，输入参数为单元编号,输出参数为包含(i,j,m)的值的一个数组i=z(1)j=z(2)m=z(3)b(i)=x(j)-x(m)c(i)=y(j)-y(m)b(j)=x(m)-x(i)c(j)=y(m)-y(i)b(m)=x(i)-x(j)c(m)=y(i)-y(j)a(i,j)=a(i,j)+1*(b(i)*b(j)+c(i)*c(j))/(4*s)!系数矩阵的计算a(i,m)=a(i,m)+1*(b(i)*b(m)+c(i)

6、*c(m))/(4*s)a(m,j)=a(m,j)+1*(b(m)*b(j)+c(m)*c(j))/(4*s)a(j,i)=a(i,j)a(j,m)=a(m,j)a(m,i)=a(i,m)a(i,i)=a(i,i)+1*(b(i)*b(i)+c(i)*c(i))/(4*s)a(m,m)=a(m,m)+1*(b(m)*b(m)+c(m)*c(m))/(4*s)a(j,j)=a(j,j)+1*(b(j)*b(j)+c(j)*c(j))/(4*s)enddor=0doi=1,6doj=7,12r(i)=r(i)-1*a(i,j)*d(j)!由于第一类边界条件，对有限元方程组的修正

7、enddoenddok=1u=50!u为求解的方程组的解，这里u=50是赋初值。w=1.22!超松弛因子count=0!count为迭代次数do!方程组的求解过程if(k==0)exitk=0doi=1,6old=u(i)new=0doj=1,6new=new+(-1)*w*a(i,j)*u(j)/a(i,i)enddonew=u(i)+new+w*r(i)/a(i,i)u(i)=newif(abs((new-old)/new)>=0.0000000001)thenk=k+1endifenddocount=coun