[Fortran] 纯文本查看 复制代码
!本程序计算可调制的正弦波激励下的氦气多个反应的介质阻挡放电
!E：电场
!f:外加正弦电压的频率
!p:气体压强
!L:求IT的时候，IT表达是前面的常数系数
!ox:时间步长
!ot:空间步长
!va:外加正弦电压
!D1，D2，DG：分别为左右两个介质层的厚度和放电间隙的宽度
!xo:空间各节点的坐标
!v0:外加正弦电压幅值
!je:电子的粒子流密度
!jp1:原子离子的粒子流密度
!jp2:分子离子的粒子流密度
!jmet:亚稳态粒子流密度
!ITz:总电流密度
!Ic:传导电流密度
!ne:电子密度
!np1:原子离子密度
!np2:分子离子密度
!De:电子的扩散系数
!Dp1:原子离子的扩散系数
!Dp2:分子离子的扩散系数
!bt1:原子离子转化系数
!bt2:分子离子转化系数
!qe:电子的电荷量
!vol:对VA求导
!miue:电子迁移率
!miup1:原子离子迁移率
!miup2:分子离子迁移率
!gama:二次电子发射系数
!s0:真空中介电常数
!sb:介质的相对介电常数
!icjifen:对传导电流在放点空间的积分
program main
implicit none
real p,L,f,ox,ot
real V,D1,D2,DG,z
real De,Dp1,Dp2
real miue,miup1,miup2
real bt1,bt2,r1,r2,et,a1,a2,a3
real qe,Tg,Te,gama,nHe
real vol,icjifen
real s0,sb,pi,kb
integer i,j,ii,k,n
real xo(10001,1)
real ne(10001,2),np1(10001,2),np2(10001,2),nmet(10001,2),nte(10001,2),E(10000,2)
real g(10000,2),R(10000,2),gg1(10000,2),RR1(10000,2),gg2(10000,2),RR2(10000,2)
real je(10000,1),jp1(10000,1),jp2(10000,1)
real XL(10000,2),XLL1(10000,2),XLL2(10000,2),XLmet(10000,2)
real IC(10000,1),ITz(1,1)
real B(9999,1),A(9998,1),C(9998,1),S(9999,1)
real BB1(9999,1),AA1(9998,1),CC1(9998,1),SS1(9999,1),BB2(9999,1),AA2(9998,1),CC2(9998,1),SS2(9999,1)
real Bmet(9999,1),Amet(9998,1),Cmet(9998,1),SSmet(9999,1)
real T(1,1),TT1(1,1),TT2(1,1),Tmet(1,1)
real vg(1,1),va(1,1), bdh1(2,1),bdh2(2,1),bdh11(2,1)
f=2.0e4
P=760.0
z=5.0e-5
ox=2.0e-5
ot=1.0e-8
L=0.22667
DG=0.2
D1=0.1
D2=0.1
V=1.6e3
De=7.2944e2
Dp1=0.5395
Dp2=0.8421
miue=4.863e2
miup1=10.4
miup2=16.7
et=2.0e-31
a1=2.31e-29
qe=1.6e-19
gama=0.01
sb=7.5
s0=8.85e-14
pi=3.14
Tg=0.02586
nHe=2.6654e19
n=0
!!!!!!!!!!!!!!!!!!!!!!空间格点的坐标
j=1
xo(1,1)=0.1
do i=2,10001
xo(i,1)=xo(i-1,1)+ox
end do
!!!!!!!!!!!!!!!!!!!!!给ne np赋初值，即t=0时刻各点ne的值
do i=1,10001
ne(i,1)=1.0e7
np1(i,1)=1.0e7
np2(i,1)=1.0e2
end do
!!!!!!!!!!!!!!!!!!!!!给E赋值，即t=0时刻各点E的值
do i=1,10000
E(i,1)=1.0e7*qe
end do
do while(.true.) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!
va(1,1)=V*sin(2.0d0*3.14d0*f*real(j)*ot)
do i=1,10000
g(i,1)=miue*ox*E(i,1)/De
R(i,1)=(exp(g(i,1))-1)/g(i,1)
gg1(i,1)=-miup1*ox*E(i,1)/Dp1
RR1(i,1)=(exp(gg1(i,1))-1)/gg1(i,1)
gg2(i,1)=-miup2*ox*E(i,1)/Dp2
RR2(i,1)=(exp(gg2(i,1))-1)/gg2(i,1)
end do
!!!!!!!!!!!!追赶法求密度
do i=1,10000
XL(i,1)=De*ot/(ox*ox*R(i,1))
XLL1(i,1)=Dp1*ot/(ox*ox*RR1(i,1))
XLL2(i,1)=Dp2*ot/(ox*ox*RR2(i,1))
end do
do ii=1,9998
B(ii,1)=1+XL(ii+1,1)+XL(ii,1)*exp(g(ii,1))
A(ii,1)=-XL(ii+1,1)
C(ii,1)=-XL(ii+1,1)*exp(g(ii+1,1))
S(ii,1)=ne(ii+1,1)+4.4*p*exp(-14.4*(p/abs(E(ii+1,1)))**0.5)*miue*abs(E(ii+1,1))*ne(ii+1,1)*ot-a1*ne(ii+1,1)*ne(ii+1,1)*np1(ii+1,1)*ot
BB1(ii,1)=1+XLL1(ii+1,1)+XLL1(ii,1)*exp(gg1(ii,1))
AA1(ii,1)=-XLL1(ii+1,1)
CC1(ii,1)=-XLL1(ii+1,1)*exp(gg1(ii+1,1))
SS1(ii,1)=np1(ii+1,1)+4.4*p*exp(-14.4*(p/abs(E(ii+1,1)))**0.5)*miue*abs(E(ii+1,1))*ne(ii+1,1)*ot-a1*ne(ii+1,1)*ne(ii+1,1)*np1(ii+1,1)*ot-et*np1(ii+1,1)*nHe*nHe*ot
BB2(ii,1)=1+XLL2(ii+1,1)+XLL2(ii,1)*exp(gg2(ii,1))
AA2(ii,1)=-XLL2(ii+1,1)
CC2(ii,1)=-XLL2(ii+1,1)*exp(gg2(ii+1,1))
SS2(ii,1)=np2(ii+1,1)+et*np1(ii+1,1)*nHe*nHe*ot
end do
B(9999,1)=1+XL(10000,1)+XL(9999,1)*exp(g(9999,1))
S(9999,1)=ne(10000,1)+4.4*p*exp(-14.4*(p/abs(E(10000,1)))**0.5)*miue*abs(E(10000,1))*ne(10000,1)*ot-a1*ne(10000,1)*ne(10000,1)*np1(10000,1)*ot
ne(1,2)=1.0e7
ne(10001,2)=1.0e7
S(1,1)=S(1,1)+XL(1,1)*ne(1,2)
S(9999,1)=S(9999,1)+XL(10000,1)*exp(g(10000,1))*ne(10001,2)!!!!!!!!!!!!!!!!!设置电子离子密度的边界值
BB1(9999,1)=1+XLL1(10000,1)+XLL1(9999,1)*exp(gg1(9999,1))
SS1(9999,1)=np1(10000,1)+4.4*p*exp(-14.4*(p/abs(E(10000,1)))**0.5)*miue*abs(E(10000,1))*ne(10000,1)*ot-a1*ne(10000,1)*ne(10000,1)*np1(10000,1)*ot
np1(1,2)=1.0d7
np1(10001,2)=1.0d7
SS1(1,1)=SS1(1,1)+XLL1(1,1)*np1(1,2)
SS1(9999,1)=SS1(9999,1)+XLL1(10000,1)*exp(gg1(10000,1))*np1(10001,2)
BB2(9999,1)=1+XLL2(10000,1)+XLL2(9999,1)*exp(gg2(9999,1))
SS2(9999,1)=np2(10000,1)+et*np1(10000,1)*nHe*nHe*ot
np2(1,2)=1.0d2
np2(10001,2)=1.0d2
SS2(1,1)=SS1(1,1)+XLL2(1,1)*np1(1,2)
!追赶法
S(1,1)=S(1,1)/B(1,1)
T=B(1,1)
SS1(1,1)=SS1(1,1)/BB1(1,1)
TT1=BB1(1,1)
SS2(1,1)=SS2(1,1)/BB2(1,1)
TT2=BB2(1,1)
do k=2,9999
B(k-1,1)=C(k-1,1)/T(1,1)
T(1,1)=B(k,1)-A(k-1,1)*B(k-1,1)
S(k,1)=(S(k,1)-A(k-1,1)*S(k-1,1))/T(1,1)
BB1(k-1,1)=CC1(k-1,1)/TT1(1,1)
TT1(1,1)=BB1(k,1)-AA1(k-1,1)*BB1(k-1,1)
SS1(k,1)=(SS1(k,1)-AA1(k-1,1)*SS1(k-1,1))/TT1(1,1)
BB2(k-1,1)=CC2(k-1,1)/TT2(1,1)
TT2(1,1)=BB2(k,1)-AA2(k-1,1)*BB2(k-1,1)
SS2(k,1)=(SS2(k,1)-AA2(k-1,1)*SS2(k-1,1))/TT2(1,1)
end do
do k=1,9998
S(9999-k,1)=S(9999-k,1)-B(9999-k,1)*S(10000-k,1)
SS1(9999-k,1)=SS1(9999-k,1)-BB1(9999-k,1)*SS1(10000-k,1)
SS2(9999-k,1)=SS2(9999-k,1)-BB2(9999-k,1)*SS2(10000-k,1)
end do
do i=2,10000
ne(i,2)=S(i-1,1)
np1(i,2)=SS1(i-1,1)
np2(i,2)=SS2(i-1,1)
end do
!过河拆桥
ne(:,1)=ne(:,2)
np1(:,1)=np1(:,2)
np2(:,1)=np2(:,2)
if(vg(1,1)>=0)then
do i=1,10000
je(i,1)=De*(ne(i,1)-ne(i+1,1)*exp(g(i,1)))/(ox*R(i,1))
end do
do i=1,10000
jp1(i,1)=Dp1*(np1(i,1)-exp(gg1(i,1))*np1(i+1,1))/(ox*RR1(i,1))
end do
do i=1,10000
jp2(i,1)=Dp2*(np2(i,1)-exp(gg2(i,1))*np2(i+1,1))/(ox*RR2(i,1))
end do
if(jp1(10000,1)+jp2(10000,1)>0)then
gama=0.01
else
gama=0.0
end if
je(10000,1)=je(10000,1)-gama*(jp1(10000,1)+jp2(10000,1))!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!je jp的边界条件
else
do i=1,10000
je(i,1)=De*(ne(i,1)-ne(i+1,1)*exp(g(i,1)))/(ox*R(i,1))
end do
do i=1,10000
jp1(i,1)=Dp1*(np1(i,1)-exp(gg1(i,1))*np1(i+1,1))/(ox*RR1(i,1))
end do
do i=1,10000
jp2(i,1)=Dp2*(np2(i,1)-exp(gg2(i,1))*np2(i+1,1))/(ox*RR2(i,1))
end do
if(jp1(10000,1)+jp2(10000,1)<0)then
gama=0.01
else
gama=0.0
end if
je(1,1)=je(1,1)-gama*(jp1(1,1)+jp2(1,1))
end if
iC=qe*(jp1+jp2-je)
icjifen=ox*sum(iC)
vol=s0*(V*sin(2.0*pi*f*real(j)*ot)-V*sin(2.0*pi*f*real(j-1)*ot))/ot
ITz(1,1)=(icjifen+vol)/L
do i=1,10000
E(i,1)=ot*(ITz(1,1)-iC(i,1))/s0+E(i,1)
end do
vg(1,1)=ox*sum(E(:,1))
bdh1(2,1)=(bdh1(1,1)+ic(1,1)*ot)
bdh2(2,1)=(bdh2(1,1)+ic(10000,1)*ot)
bdh1(1,1)=bdh1(2,1)
bdh11(1,1)=-bdh1(1,1)
bdh2(1,1)=bdh2(2,1)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
if(mod(j,100)==0)then
open(1,file="ITz.dat")
open(2,file="vol.dat")
open(3,file="vg.dat")
open(4,file="bdh11.dat")
open(5,file="bdh2.dat")
open(6,file="ic.dat")
write(1,*) j,ITz(1,1)
write(2,*) j,va(1,1)
write(3,*) j,vg(1,1)
write(4,*) j,bdh11(1,1)
write(5,*) j,bdh2(1,1)
write(6,*) j,ic(1,1)
end if
j=j+1
write(*,*)j
if(mod(j,5000)==0)then
n=n+1
open(10,file="n.dat")
write(10,*) n
end if
end do
close(1)
close(2)
close(3)
close(4)
close(5)
close(6)
close(10)
stop
end program