计算电流线圈产生的磁场1.doc

（二）求截面为矩形的圆线圈周围产生的磁场一、数值方法 （一）数学模型：所研究的电流圆线圈产生磁场的问题在柱坐标系下研究， 根据磁场强度跟矢势之间的关系，得到磁场； 磁场为，矢势为（由具有轴对称得到）所以在柱坐标系中，由公式-得即，（1）先求矢势一个电流为，半径为的线圆环周围空间产生的磁场，其矢势表示为推广到截面为矩形的圆环线圈中其中为矩形截面的面积，为矩形截面的两边距圆环中心的距离，为矩形截面的上下面的轴坐标。 （二）数值模型离散化（均匀网格有限差分）（1） 高斯方法计算三重积分（参考书:徐士良常用算法程序集第二版）（2） 根据一阶三点公式来求解磁场的分布： 磁场分量,对于点处二 数值计算程序SUBROUTINE FGAUS(N,JS,X,FS,F,S,Z,R)DIMENSION JS(N),X(N)DIMENSION T(5),C(5),D(2,11),CC(11),IS(2,11)DATA T/-0.90617,-0.538469,0.0,0.538469,0.90617/DATA C/0.2369,0.47862,0.568889,0.4786,0.2369/M=1D(1,N+1)=1.0D(2,N+1)=1.010 DO 20 J=M,NCALL FS (J,N,X,DN,UP)D(1,J)=0.5*(UP-DN)/JS(J)CC(J)=D(1,J)+DNX(J)=D(1,J)*T(1)+CC(J)D(2,J)=0.0IS(1,J)=1IS(2,J)=120 CONTINUEJ=N30 K=IS(1,J)IF(J.EQ.N)THENP=F(N,X,Z,R)ELSEP=1.0ENDIFD(2,J)=D(2,J+1)*D(1,J+1)*P*C(K)+D(2,J)IS(1,J)=IS(1,J)+1IF(IS(1,J).GT.5)THEN IF(IS(2,J).GE.JS(J) THEN J=J-1 IF(J.EQ.0) THEN S=D(2,1)*D(1,1) RETURN ENDIF GOTO 30 ENDIF IS(2,J)=IS(2,J)+1 CC(J)=CC(J)+D(1,J)*2.0 IS(1,J)=1 ENDIF K=IS(1,J) X(J)=D(1,J)*T(K)+CC(J) IF(J.EQ.N) GOTO 30 M=J+1 GOTO 10 ENDEXTERNAL FS,FDIMENSION JS(3),X(3),Z(30),R(30),a0(30,30),BB(30,30),BR(30,30),BZ(30,30)DATA JS/4,4,4/N=3c=1H=0.5Z(1)=1.5DO I=1,29Z(I+1)=Z(I)+HENDDOR(1)=1.5do j=1,29R(J+1)=R(J)+HENDDODO I=1,30 DO J=1,30CALL FGAUS(N,JS,X,FS,F,S,Z(I),R(J)a0(i,j)=c*sWRITE(1,*)Z(I),R(J),SOPEN(1,FILE=DUHAI.DAT)ENDDOENDDO!求解磁场的Brdo i=1,30Br(i,1)=(-3*A0(i,1)+4*A0(i,2)-A0(i,3)/2*hBr(i,30)=(A0(i,28)-4*A0(i,29)+3*A0(i,30)/2*hdo j=2,29Br(i,j)=(A0(i,j-1)-A0(i,j+1)/2*hend doend do!求解磁场的Bzdo j=1,30Bz(1,j)=(-3*A0(1,j)+4*A0(2,j)-A0(3,j)/2*h+A0(1,j)/r(1)Bz(30,j)=(A0(28,j)-4*A0(29,j)+3*A0(30,j)/2*h+A0(30,j)/r(30)do i=2,29Bz(i,j)=(A0(i+1,j)-A0(i-1,j)/2*h+A0(i,j)/r(i)end doend dodo i=1,30 do j=1,30BB(i,j)=sqrt(Br(i,j)*2+Bz(i,j)*2)end doend dodo i=1,30 do j=1,30 write(*,*)z(i),r(j), bb(i,j) write(2,*)z(i),r(j), bb(i,j)open (2,file=cichang.dat) enddoEnddoEndSUBROUTINE FS(J,N,X,DN,UP)DIMENSION X(N)IF(J.EQ.1) THENDN=0.5UP=1.0ELSEIF(J.EQ.2)THENDN=0.5UP=1.0ELSEIF(J.EQ.3) THENDN=0.0 UP=2*3.1415926 ENDIF RETURN END FUNCTION F(N,X,Z,R) DIMENSION X(N) f=x(2)*cos(x(1)/sqrt(z-x(3)*2+R*2+x(2)*2-2*R*x