JGonzalez/vlaplex
1
0
Fork 0
Code Issues Pull requests Projects Releases Packages Wiki Activity Actions

Very small progress, E at r0 still high

Browse source
This commit is contained in:
Jorge Gonzalez 2024-11-12 11:37:40 +01:00
commit a0b73194d2
1 changed files with 10 additions and 8 deletions
Download patch file Download diff file Expand all files Collapse all files

18
plasmaExpansion.f90
View file

@ -165,8 +165,8 @@ program plasmaExpansion
u_i = 0.0_dp u_i = 0.0_dp
E_i = 0.0_dp E_i = 0.0_dp
T_i = 0.0_dp T_i = 0.0_dp
n_e = 1.0_dp n_e = 0.0_dp
T_e = 1.0_dp T_e = 0.0_dp
Zave = 0.0_dp Zave = 0.0_dp
phi = 0.0_dp phi = 0.0_dp
phi_old = 0.0_dp phi_old = 0.0_dp
@ -209,6 +209,7 @@ program plasmaExpansion
! Set boundary values ! Set boundary values
phi0 = 0.0_dp / phi_ref ! Dirichlet phi0 = 0.0_dp / phi_ref ! Dirichlet
phi(1) = phi0 ! Dirichlet
! phi0 = phi(1) ! Neumann ! phi0 = phi(1) ! Neumann
! phiF = 0.0_dp ! Dirichlet ! phiF = 0.0_dp ! Dirichlet
allocate(f0(j0:nv)) allocate(f0(j0:nv))
@ -232,6 +233,7 @@ program plasmaExpansion
! f0(j0:nv) = v(j0:nv)**2 / sqrt(PI*Temp_bc**3) * exp(-(v(j0:nv) - u_bc)**2 / Temp_bc) ! f0(j0:nv) = v(j0:nv)**2 / sqrt(PI*Temp_bc**3) * exp(-(v(j0:nv) - u_bc)**2 / Temp_bc)
f0(j0:nv) = 1.0_dp / sqrt(PI*Temp_bc) * exp(-(v(j0:nv) - u_bc)**2 / Temp_bc) f0(j0:nv) = 1.0_dp / sqrt(PI*Temp_bc) * exp(-(v(j0:nv) - u_bc)**2 / Temp_bc)
f0 = f0 * n_bc / (sum(f0)*dv) f0 = f0 * n_bc / (sum(f0)*dv)
T_e = Temp_bc
! Boundary conditions ! Boundary conditions
! r = r0, v>0 ! r = r0, v>0
@ -280,6 +282,11 @@ program plasmaExpansion
! Store previous value ! Store previous value
phi_old = phi phi_old = phi
! Calculate distribution of electrons
! n_e = Zave * n_i ! Quasi-neutral
n_e = Zave(1) * n_i(1) * exp((phi_old - phi0) / T_e) ! Isothermal (Boltzmann)
! n_e = Zave(1) * n_i(1) * (1.0_dp + ((gamma_e - 1.0_dp)/gamma_e*(phi_old-phi0)/T_e)**(1.0_dp/(gamma_e - 1.0_dp))) ! Not working
! Diagonal matrix for Newton integration scheme ! Diagonal matrix for Newton integration scheme
db_dphi = n_e / T_e ! Isotropic db_dphi = n_e / T_e ! Isotropic
diag = -2.0_dp / dr**2 - db_dphi diag = -2.0_dp / dr**2 - db_dphi
@ -306,12 +313,7 @@ program plasmaExpansion
! Iterate system ! Iterate system
call dgtsv(nr, 1, diag_low, diag, diag_high, Res, nr, info) call dgtsv(nr, 1, diag_low, diag, diag_high, Res, nr, info)
phi = phi_old + 1.0e-1_dp*Res phi = phi_old + 1.0e-1_dp*Res
phi0=phi(1) ! Neumann ! phi0=phi(1) ! Neumann
! Calculate distribution of electrons
T_e = Temp_bc
n_e = Zave(1) * n_i(1) * exp((phi - phi0) / T_e)
! n_e = Zave(1) * n_i(1) * (1.0_dp + ((gamma_e - 1.0_dp)/gamma_e*(phi-phi0)/T_e)**(1.0_dp/(gamma_e - 1.0_dp))) ! Not working
! Check if the solution has converged ! Check if the solution has converged
phiConv = maxval(abs(Res),1) phiConv = maxval(abs(Res),1)