[Fortran] 纯文本查看 复制代码
program main
implicit none
real p,V0,f
!real RHS
real ox,ot
real De,Dp,qe,bt,KB
!real RS,git,rit
real miue,miup
real icjifen,vol,gama
real s0, sb,d1,d2,l,Dg
integer i,j,ii,k,error
!real am,bm,cm,dm,em
real xo(1001,1)
real ne(1001,2),np(1001,2),E(1000,2)
real g(1000,2),R(1000,2),gg(1000,2),RR(1000,2) ,gte(1000,2),Rte(1000,2)
real je(1000,1),jp(1000,1),jq(1000,1),XL(1000,2),XLL(1000,2),Xte(1000,2)
real IC(1000,1)
real B(999,1),A(998,1),C(998,1),S(999,1)
real Bte(999,1),Ate(998,1),Cte(998,1),Ste(999,1)
real BB(999,1),AA(998,1),CC(998,1),SS(999,1),T(1,1),TT(1,1),ITz(1,1) ,Tte(1,1)
real vg(1,1),va(1,1),bdh1(2,1),bdh2(2,1),bdh11(2,1)
real Te(1000,1)
V0=1.5e3
P=760.0
F=3.3e4
s0=8.85E-14
bt=1.389e-28
sb=7.5
D1=0.1
D2=0.1
Dg=0.2
L=0.2267
KB=1.380605e-23
gama=0.01
ox=2.0E-4
ot=1.0E-9
miue=4.863e2
De=7.2944e2
miup=10.4
Dp=0.53947
qe=1.6E-19
!!!!!!!!!!!!!!!!!!!!!!空间格点的坐标
j=1
xo(1,1)=0.1d0
do i=2,1001
xo(i,1)=xo(i-1,1)+ox
end do
!!!!!!!!!!!!!!!!!!!!!给ne np赋初值，即t=0时刻各点ne的值
do i=1,1001
ne(i,1)=1.0d7
np(i,1)=1.0d7
te(i,1)=0.026
end do
!!!!!!!!!!!!!!!!!!!!!给E赋值，即t=0时刻各点E的值
do i=1,1000
E(i,1)= 1.0d7*qe
end do
do while(.true.) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!
va(1,1)=V0*sin(2.0d0*3.14d0*f*real(j)*ot)
do i=1,1000
g(i,1)=miue*ox*E(i,1)/De
R(i,1)=(exp(g(i,1))-1)/g(i,1)
gg(i,1)=-miup*ox*E(i,1)/Dp
RR(i,1)=(exp(gg(i,1))-1)/gg(i,1)
gTE(i,1)=miue*ox*E(i,1)/De
RTE(i,1)=(exp(gTE(i,1))-1)/gTE(i,1)
end do
!!!!!!!!!!!!追赶法求密度
do i=1,1000
XL(i,1)=De*ot/(ox*ox*R(i,1))
XLL(i,1)=Dp*ot/(ox*ox*RR(i,1))
XTE(i,1)=5*De*ot /(3*ox*ox*RTE(i,1))
end do
do ii=1,998
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.0*(p/abs(E(ii+1,1)))**0.5)*miue*abs(E(ii+1,1))*ot*ne(ii+1,1)-bt*ot*ne(ii+1,1)*ne(ii+1,1)*np(ii+1,1)
BB(ii,1)=1+XLL(ii+1,1)+XLL(ii,1)*exp(gg(ii,1))
AA(ii,1)=-XLL(ii+1,1)
CC(ii,1)=-XLL(ii+1,1)*exp(gg(ii+1,1))
SS(ii,1)=np(ii+1,1)+4.4*p*exp(-14.0*(p/abs(E(ii+1,1)))**0.5)*miue*abs(E(ii+1,1))*ot*ne(ii+1,1)-bt*ot*ne(ii+1,1)*ne(ii+1,1)*np(ii+1,1)
BTE(ii,1)=1+XTE(ii+1,1)+XTE(ii,1)*exp(gTE(ii,1))
ATE(ii,1)=-XTE(ii+1,1)
CTE(ii,1)=-XTE(ii+1,1)*exp(gTE(ii+1,1))
STE(ii,1)=ne(ii+1,1)-qe*abs(E(ii+1,1))*je(ii+1,1)*ot-bt*ot*ne(ii+1,1)*ne(ii+1,1)*np(ii+1,1)*(3*KB*TE(ii+1,1)/2 )
end do
end do
B(999,1)=1+XL(1000,1)+XL(999,1)*exp(g(999,1))
S(999,1)=ne(1000,1)+4.4*p*exp(-14.0*(p/abs(E(1000,1)))**0.5)*miue*abs(E(1000,1))*ot*ne(1000,1)-bt*ot*ne(1000,1)*ne(1000,1)*np(1000,1)
ne(1,2)=1.0d7
ne(1001,2)=1.0d7
S(1,1)=S(1,1)+XL(1,1)*ne(1,2)
S(999,1)=S(999,1)+XL(1000,1)*exp(g(1000,1))*ne(1001,2)!!!!!!!!!!!!!!!!!设置电子离子密度的边界值
BB(999,1)=1+XLL(1000,1)+XLL(999,1)*exp(gg(999,1))
SS(999,1)=np(1000,1)+4.4*p*exp(-14.0*(p/abs(E(1000,1)))**0.5)*miue*abs(E(1000,1))*ot*ne(1000,1)-bt*ot*ne(1000,1)*ne(1000,1)*np(1000,1)
np(1,2)=np(2,1)
np(1001,2)=np(1000,1)
SS(1,1)=SS(1,1)+XLL(1,1)*np(1,2)
SS(999,1)=SS(999,1)+XLL(1000,1)*exp(gg(1000,1))*np(1001,2)
BTE(999,1)=1+XTE(1000,1)+XTE(999,1)*exp(gTE(999,1))
STE(999,1)=-qe*abs(E(1000,1))*ot*je(1000,1)-bt*ot*ne(1000,1)*ne(1000,1)*np(1000,1)*3*KB*TE(1000,1)/2
TE(1,2)=1.5
TE(1001,2)=1.5
STE(1,1)=STE(1,1)+XTE(1,1)*TE(1,2)
STE(999,1)=STE(999,1)+XTE(1000,1)*exp(gTE(1000,1))*TE(1001,2)!!!!!!!!
!追赶法
S(1,1)=S(1,1)/B(1,1)
T=B(1,1)
SS(1,1)=SS(1,1)/BB(1,1)
TT=BB(1,1)
STE(1,1)=STE(1,1)/BTE(1,1)
TTE=BTE(1,1)
do k=2,999
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)
BB(k-1,1)=CC(k-1,1)/TT(1,1)
TT(1,1)=BB(k,1)-AA(k-1,1)*BB(k-1,1)
SS(k,1)=(SS(k,1)-AA(k-1,1)*SS(k-1,1))/TT(1,1)
BTE(k-1,1)=CTE(k-1,1)/TTE(1,1)
TTE(1,1)=BTE(k,1)-ATE(k-1,1)*BTE(k-1,1)
STE(k,1)=(STE(k,1)-ATE(k-1,1)*STE(k-1,1))/TTE(1,1)
end do
do k=1,998
S(999-k,1)=S(999-k,1)-B(999-k,1)*S(1000-k,1)
SS(999-k,1)=SS(999-k,1)-BB(999-k,1)*SS(1000-k,1)
STE(999-k,1)=STE(999-k,1)-BTE(999-k,1)*STE(1000-k,1)
end do
do i=2,1000
ne(i,2)=S(i-1,1)
np(i,2)=SS(i-1,1)
TE(i,2)=STE(i-1,1)
end do
!过河拆桥
ne(:,1)=ne(:,2)
np(:,1)=np(:,2)
TE(:,1)=TE(:,2)
if(vg(1,1)>=0)then
do i=1,1000
je(i,1)=De*(ne(i,1)-ne(i+1,1)*exp(g(i,1)))/(ox*R(i,1))
jp(i,1)=Dp*(np(i,1)-exp(gg(i,1))*np(i+1,1))/(ox*RR(i,1))
jq(1,1)=3*(De*(ne(i,1)-ne(i+1,1)*exp(gte(i,1))))/5*(ox*Rte(i,1))
end do
if(jp(1000,1)>0)then
gama=0.01d0
else
gama=0.0d0
end if
je(1000,1)=je(1000,1)-gama*jp(1000,1)!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!je jp的边界条件
else
do i=1,1000
je(i,1)=De*(ne(i,1)-ne(i+1,1)*exp(g(i,1)))/(ox*R(i,1))
jq(i,1)= 3*(De*(ne(i,1)-ne(i+1,1)*exp(gte(i,1))))/5*(ox*Rte(i,1))
jp(i,1)=Dp*(np(i,1)-exp(gg(i,1))*np(i+1,1))/(ox*RR(i,1))
end do
if(jp(1000,1)<0)then
gama=0.01d0
else
gama=0.0d0
end if
je(1,1)=je(1,1)-gama*jp(1,1)
end if
iC=qe*(jp-je)
icjifen=ox*sum(iC)
vol=s0*(V0*sin(2.0d0*3.14d0*f*real(j)*ot)-V0*sin(2.0d0*3.14d0*f*real(j-1)*ot))/ot
ITz(1,1)=(vol+icjifen)/L
do i=1,1000
E(i,1)=ot*(ITz(1,1)-iC(i,1))/s0+E(i,1)
end do
vg(1,1)=ox*sum(E(:,1))!!!!!2时刻的气体电压
bdh1(2,1)=(bdh1(1,1)+ic(1,1)*ot)
bdh2(2,1)=(bdh2(1,1)+ic(1000,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="va.dat")
open(3,file="vg.dat")
open(4,file="bdh11.dat")
open(5,file="bdh2.dat")
open(6,file="Te.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,Te(1,1)
end if
j=j+1
write(*,*)j
close(1)
close(2)
close(3)
close(4)
close(5)
close(6)
stop
end program
F 币
贡献
权杖
窥视卡
发表于 2017-10-25 10:34:04
微信
收藏
点赞
点踩
提升卡
置顶卡
变色卡
千斤顶
显身卡
楼主
发表于 2017-10-26 13:12:39