求截面为矩形的圆线圈周围产生的磁场
一、数值方法
（一）数学模型：所研究的电流圆线圈产生磁场的问题在柱坐标系下研究， 根据磁场强度跟矢势之间的关系，得到磁场；
磁场为B ，矢势为A
B A =∇⨯
r r z z A A e A e A e θθ=++ A e θθ= (,)A r z e θθ=
（由A 具有轴对称得到）
所以B A =∇⨯ A e θθ=∇⨯
在柱坐标系中，由公式1()()11()()r r z z
z r
r z r
z f f e f e f e f f f r z
f f f z r
f
f rf r r r θθθ
θθθθ
⎧=++⎪∂∂⎪∇⨯=-⎪∂∂⎪⎨∂∂∇⨯=-⎪∂∂⎪∂∂
⎪∇⨯=-⎪∂∂⎩
-得
B A =∇⨯ 1()r z f e rf e z r r θθ∂∂
=-+∂∂
即r A B z θ
∂=-∂，1()z B rA r r θ∂
=∂
（1）先求矢势A
4L Idl A r
μπ=⎰
一个电流为I ，半径为a 的线圆环周围空间产生的磁场，其矢势表示为
202220cos (,)42cos Ia A r z d r z a ar πθμϕ
ϕπϕ=++-⎰
推广到截面为矩形的圆环线圈中
22
11202220
cos (,)4()2cos R z R z I r A r z d dz dr s r z z r r r πθμϕ
ϕπϕ'''='''+-+-⎰⎰⎰
其中S 为矩形截面的面积，12,R R 为矩形截面的两边距圆环中心的距离，12,z z 为矩形截面的上下面的z 轴坐标。

（二）数值模型离散化（均匀网格有限差分）
（1）高斯方法计算三重积分（参考书:徐士良常用算法程序集第二版）
（2）根据一阶三点公式来求解磁场的分布：
磁场分量r B ,z B
,1,1()()2i j i j r A A A B z z
θθθ-+-∂=-=∂∆ 1,1,,,()()()1()2i j i j i j z i j A A A B rA r r r r θθθθ+--∂=
=+∂∆ 对于0r →点处
1,1,0
00()()1lim lim ()lim()2i j i j z r r r A A A A A B rA r r r r r r θθθθθθ+-→→→-∂∂∂==+==∂∂∂∆ 22z r B B B +=
二 数值计算程序
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=1
D(1,N+1)=1.0
D(2,N+1)=1.0
10 DO 20 J=M,N
CALL FS (J,N,X,DN,UP)
D(1,J)=0.5*(UP-DN)/JS(J)
CC(J)=D(1,J)+DN
X(J)=D(1,J)*T(1)+CC(J)
D(2,J)=0.0
IS(1,J)=1
IS(2,J)=1
20 CONTINUE
J=N
30 K=IS(1,J)
IF(J.EQ.N)THEN
P=F(N,X,Z,R)
ELSE
P=1.0
ENDIF
D(2,J)=D(2,J+1)*D(1,J+1)*P*C(K)+D(2,J) IS(1,J)=IS(1,J)+1
IF(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
END
EXTERNAL FS,F
DIMENSION
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=3
c=1
H=0.5
Z(1)=1.5
DO I=1,29
Z(I+1)=Z(I)+H
ENDDO
R(1)=1.5
do j=1,29
R(J+1)=R(J)+H
ENDDO
DO I=1,30
DO J=1,30
CALL FGAUS(N,JS,X,FS,F,S,Z(I),R(J))
a0(i,j)=c*s
WRITE(1,*)Z(I),R(J),S
OPEN(1,FILE='DUHAI.DAT')
ENDDO
ENDDO
!求解磁场的Br
do i=1,30
Br(i,1)=(-3*A0(i,1)+4*A0(i,2)-A0(i,3))/2*h
Br(i,30)=(A0(i,28)-4*A0(i,29)+3*A0(i,30))/2*h
do j=2,29
Br(i,j)=(A0(i,j-1)-A0(i,j+1))/2*h
end do
end do
!求解磁场的Bz
do j=1,30
Bz(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,29
Bz(i,j)=(A0(i+1,j)-A0(i-1,j))/2*h+A0(i,j)/r(i)
end do
end do
do i=1,30
do j=1,30
BB(i,j)=sqrt(Br(i,j)**2+Bz(i,j)**2)
end do
end do
do 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')
enddo
Enddo
End
SUBROUTINE FS(J,N,X,DN,UP)
DIMENSION X(N)
IF(J.EQ.1) THEN
DN=0.5
UP=1.0
ELSEIF(J.EQ.2)THEN
DN=0.5
UP=1.0
ELSEIF(J.EQ.3) THEN
DN=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(2)*cos(x(1))) RETURN
END
三计算结果矢势A的分布如下：
磁场B的分布如下
四结果讨论
由上面的计算的结果可以看出磁场是随着z
r,的增加而减小，反之，随着z
r,减
小而增加，当增加到一定的程度上，几乎是趋向于零；。