[其他行业算法] 如何将电子能量守恒方程加进去，求出电子温度

!本程序计算可调制的正弦波激励下的氦气多个反应的介质阻挡放电
!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