如何将电子能量守恒方程加进去，求出电子温度

Nimo

[Fortran] 纯文本查看 复制代码

查看代码

打印?

001	!本程序计算可调制的正弦波激励下的氦气多个反应的介质阻挡放电

002

!E：电场

003	!f:外加正弦电压的频率

004

!p:气体压强

005	!L:求IT的时候，IT表达是前面的常数系数

006

!ox:时间步长

007

!ot:空间步长

008	!va:外加正弦电压

009	!D1，D2，DG：分别为左右两个介质层的厚度和放电间隙的宽度

010	!xo:空间各节点的坐标

011	!v0:外加正弦电压幅值

012	!je:电子的粒子流密度

013	!jp1:原子离子的粒子流密度

014	!jp2:分子离子的粒子流密度

015	!jmet:亚稳态粒子流密度

016	!ITz:总电流密度

017	!Ic:传导电流密度

018

!ne:电子密度

019	!np1:原子离子密度

020	!np2:分子离子密度

021	!De:电子的扩散系数

022	!Dp1:原子离子的扩散系数

023	!Dp2:分子离子的扩散系数

024	!bt1:原子离子转化系数

025	!bt2:分子离子转化系数

026	!qe:电子的电荷量

027	!vol:对VA求导

028	!miue:电子迁移率

029	!miup1:原子离子迁移率

030	!miup2:分子离子迁移率

031	!gama:二次电子发射系数

032	!s0:真空中介电常数

033	!sb:介质的相对介电常数

034	!icjifen:对传导电流在放点空间的积分

035

036

037	program main

038	implicit none

039	real p,L,f,ox,ot

040	real V,D1,D2,DG,z

041	real De,Dp1,Dp2

042	real miue,miup1,miup2

043	real bt1,bt2,r1,r2,et,a1,a2,a3

044	real qe,Tg,Te,gama,nHe

045	real vol,icjifen

046	real s0,sb,pi,kb

047	integer i,j,ii,k,n

048

049	real xo(10001,1)

050	real ne(10001,2),np1(10001,2),np2(10001,2),nmet(10001,2),nte(10001,2),E(10000,2)

051	real g(10000,2),R(10000,2),gg1(10000,2),RR1(10000,2),gg2(10000,2),RR2(10000,2)

052	real je(10000,1),jp1(10000,1),jp2(10000,1)

053	real XL(10000,2),XLL1(10000,2),XLL2(10000,2),XLmet(10000,2)

054	real IC(10000,1),ITz(1,1)

055	real B(9999,1),A(9998,1),C(9998,1),S(9999,1)

056	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)

057	real Bmet(9999,1),Amet(9998,1),Cmet(9998,1),SSmet(9999,1)

058	real T(1,1),TT1(1,1),TT2(1,1),Tmet(1,1)

059	real vg(1,1),va(1,1), bdh1(2,1),bdh2(2,1),bdh11(2,1)

060

061

062

f=2.0e4

063

P=760.0

064

z=5.0e-5

065

ox=2.0e-5

066

ot=1.0e-8

067

L=0.22667

068

DG=0.2

069

D1=0.1

070

D2=0.1

071

V=1.6e3

072

073	De=7.2944e2

074	Dp1=0.5395

075	Dp2=0.8421

076	miue=4.863e2

077	miup1=10.4

078	miup2=16.7

079

080	et=2.0e-31

081	a1=2.31e-29

082

083	qe=1.6e-19

084

gama=0.01

085

sb=7.5

086	s0=8.85e-14

087

pi=3.14

088	Tg=0.02586

089	nHe=2.6654e19

090

n=0

091

092	!!!!!!!!!!!!!!!!!!!!!!空间格点的坐标

093

j=1

094	xo(1,1)=0.1

095	do i=2,10001

096	xo(i,1)=xo(i-1,1)+ox

097

end do

098	!!!!!!!!!!!!!!!!!!!!!给ne np赋初值，即t=0时刻各点ne的值

099	do i=1,10001

100	ne(i,1)=1.0e7

101	np1(i,1)=1.0e7

102	np2(i,1)=1.0e2

103

end do

104

105	!!!!!!!!!!!!!!!!!!!!!给E赋值，即t=0时刻各点E的值

106	do i=1,10000

107	E(i,1)=1.0e7*qe

108

end do

109

110	do while(.true.)   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

111	!!!!!!!!!!!!!!!!!!!!

112

113	va(1,1)=V*sin(2.0d0*3.14d0*f*real(j)*ot)

114

115

116

117

118	do i=1,10000

119	g(i,1)=miue*ox*E(i,1)/De

120	R(i,1)=(exp(g(i,1))-1)/g(i,1)

121	gg1(i,1)=-miup1*ox*E(i,1)/Dp1

122	RR1(i,1)=(exp(gg1(i,1))-1)/gg1(i,1)

123	gg2(i,1)=-miup2*ox*E(i,1)/Dp2

124	RR2(i,1)=(exp(gg2(i,1))-1)/gg2(i,1)

125

126

end do

127

128

129

130	!!!!!!!!!!!!追赶法求密度

131	do i=1,10000

132	XL(i,1)=De*ot/(ox*ox*R(i,1))

133	XLL1(i,1)=Dp1*ot/(ox*ox*RR1(i,1))

134	XLL2(i,1)=Dp2*ot/(ox*ox*RR2(i,1))

135

end do

136

137

138	do ii=1,9998

139	B(ii,1)=1+XL(ii+1,1)+XL(ii,1)*exp(g(ii,1))

140	A(ii,1)=-XL(ii+1,1)

141	C(ii,1)=-XL(ii+1,1)*exp(g(ii+1,1))

142	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

143

144	BB1(ii,1)=1+XLL1(ii+1,1)+XLL1(ii,1)*exp(gg1(ii,1))

145	AA1(ii,1)=-XLL1(ii+1,1)

146	CC1(ii,1)=-XLL1(ii+1,1)*exp(gg1(ii+1,1))

147	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

148

149	BB2(ii,1)=1+XLL2(ii+1,1)+XLL2(ii,1)*exp(gg2(ii,1))

150	AA2(ii,1)=-XLL2(ii+1,1)

151	CC2(ii,1)=-XLL2(ii+1,1)*exp(gg2(ii+1,1))

152	SS2(ii,1)=np2(ii+1,1)+et*np1(ii+1,1)*nHe*nHe*ot

153

154

155

156

end do

157

158

159	B(9999,1)=1+XL(10000,1)+XL(9999,1)*exp(g(9999,1))

160	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

161	ne(1,2)=1.0e7

162	ne(10001,2)=1.0e7

163	S(1,1)=S(1,1)+XL(1,1)*ne(1,2)

164	S(9999,1)=S(9999,1)+XL(10000,1)*exp(g(10000,1))*ne(10001,2)!!!!!!!!!!!!!!!!!设置电子离子密度的边界值

165

166	BB1(9999,1)=1+XLL1(10000,1)+XLL1(9999,1)*exp(gg1(9999,1))

167	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

168	np1(1,2)=1.0d7

169	np1(10001,2)=1.0d7

170	SS1(1,1)=SS1(1,1)+XLL1(1,1)*np1(1,2)

171	SS1(9999,1)=SS1(9999,1)+XLL1(10000,1)*exp(gg1(10000,1))*np1(10001,2)

172

173	BB2(9999,1)=1+XLL2(10000,1)+XLL2(9999,1)*exp(gg2(9999,1))

174	SS2(9999,1)=np2(10000,1)+et*np1(10000,1)*nHe*nHe*ot

175	np2(1,2)=1.0d2

176	np2(10001,2)=1.0d2

177	SS2(1,1)=SS1(1,1)+XLL2(1,1)*np1(1,2)

178

179

180

!追赶法

181	S(1,1)=S(1,1)/B(1,1)

182

T=B(1,1)

183

184	SS1(1,1)=SS1(1,1)/BB1(1,1)

185	TT1=BB1(1,1)

186

187	SS2(1,1)=SS2(1,1)/BB2(1,1)

188	TT2=BB2(1,1)

189

190

191	do k=2,9999

192

193	B(k-1,1)=C(k-1,1)/T(1,1)

194	T(1,1)=B(k,1)-A(k-1,1)*B(k-1,1)

195	S(k,1)=(S(k,1)-A(k-1,1)*S(k-1,1))/T(1,1)

196

197	BB1(k-1,1)=CC1(k-1,1)/TT1(1,1)

198	TT1(1,1)=BB1(k,1)-AA1(k-1,1)*BB1(k-1,1)

199	SS1(k,1)=(SS1(k,1)-AA1(k-1,1)*SS1(k-1,1))/TT1(1,1)

200

201	BB2(k-1,1)=CC2(k-1,1)/TT2(1,1)

202	TT2(1,1)=BB2(k,1)-AA2(k-1,1)*BB2(k-1,1)

203	SS2(k,1)=(SS2(k,1)-AA2(k-1,1)*SS2(k-1,1))/TT2(1,1)

204

205

end do

206

207	do k=1,9998

208	S(9999-k,1)=S(9999-k,1)-B(9999-k,1)*S(10000-k,1)

209	SS1(9999-k,1)=SS1(9999-k,1)-BB1(9999-k,1)*SS1(10000-k,1)

210	SS2(9999-k,1)=SS2(9999-k,1)-BB2(9999-k,1)*SS2(10000-k,1)

211

end do

212

213	do i=2,10000

214	ne(i,2)=S(i-1,1)

215	np1(i,2)=SS1(i-1,1)

216	np2(i,2)=SS2(i-1,1)

217

end do

218

219

220

!过河拆桥

221	ne(:,1)=ne(:,2)

222	np1(:,1)=np1(:,2)

223	np2(:,1)=np2(:,2)

224

225

226

227	if(vg(1,1)>=0)then

228	do i=1,10000

229	je(i,1)=De*(ne(i,1)-ne(i+1,1)*exp(g(i,1)))/(ox*R(i,1))

230

end do

231	do i=1,10000

232	jp1(i,1)=Dp1*(np1(i,1)-exp(gg1(i,1))*np1(i+1,1))/(ox*RR1(i,1))

233

end do

234	do i=1,10000

235	jp2(i,1)=Dp2*(np2(i,1)-exp(gg2(i,1))*np2(i+1,1))/(ox*RR2(i,1))

236

end do

237

238	if(jp1(10000,1)+jp2(10000,1)>0)then

239

gama=0.01

240

else

241

gama=0.0

242

end if

243	je(10000,1)=je(10000,1)-gama*(jp1(10000,1)+jp2(10000,1))!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!je jp的边界条件

244

245

else

246	do i=1,10000

247	je(i,1)=De*(ne(i,1)-ne(i+1,1)*exp(g(i,1)))/(ox*R(i,1))

248

end do

249	do i=1,10000

250	jp1(i,1)=Dp1*(np1(i,1)-exp(gg1(i,1))*np1(i+1,1))/(ox*RR1(i,1))

251

end do

252

253	do i=1,10000

254	jp2(i,1)=Dp2*(np2(i,1)-exp(gg2(i,1))*np2(i+1,1))/(ox*RR2(i,1))

255

end do

256

257	if(jp1(10000,1)+jp2(10000,1)<0)then

258

gama=0.01

259

else

260

gama=0.0

261

end if

262	je(1,1)=je(1,1)-gama*(jp1(1,1)+jp2(1,1))

263

264

end if

265

266	iC=qe*(jp1+jp2-je)

267	icjifen=ox*sum(iC)

268

269	vol=s0*(V*sin(2.0*pi*f*real(j)*ot)-V*sin(2.0*pi*f*real(j-1)*ot))/ot

270

271

272

273

274	ITz(1,1)=(icjifen+vol)/L

275	do i=1,10000

276	E(i,1)=ot*(ITz(1,1)-iC(i,1))/s0+E(i,1)

277

end do

278

279	vg(1,1)=ox*sum(E(:,1))

280

281	bdh1(2,1)=(bdh1(1,1)+ic(1,1)*ot)

282

283	bdh2(2,1)=(bdh2(1,1)+ic(10000,1)*ot)

284

285

286	bdh1(1,1)=bdh1(2,1)

287	bdh11(1,1)=-bdh1(1,1)

288	bdh2(1,1)=bdh2(2,1)

289	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

290	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

291	if(mod(j,100)==0)then

292	open(1,file="ITz.dat")

293	open(2,file="vol.dat")

294	open(3,file="vg.dat")

295	open(4,file="bdh11.dat")

296	open(5,file="bdh2.dat")

297	open(6,file="ic.dat")

298

299	write(1,*)   j,ITz(1,1)

300	write(2,*)   j,va(1,1)

301	write(3,*)   j,vg(1,1)

302	write(4,*)   j,bdh11(1,1)

303	write(5,*)   j,bdh2(1,1)

304	write(6,*)   j,ic(1,1)

305

end if

306

307

308

309

310

311

312

313

j=j+1

314	write(*,*)j

315	if(mod(j,5000)==0)then

316

n=n+1

317	open(10,file="n.dat")

318	write(10,*) n

319

end if

320

321

end do

322

close(1)

323

close(2)

324

close(3)

325

close(4)

326

close(5)

327

close(6)

328

329

close(10)

330

stop

331	end program